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Preface 


This document contains the proceedings of the Training Workshop on Multiscale 
Modeling, Simulation and Visualization and Their Potential for Future Aerospace 
Systems held at NASA Langley Research Center, Hampton, Virginia, March 5 -6, 2002. 
The workshop was jointly sponsored by Old Dominion University and NASA. 
Workshop attendees came from NASA, other government agencies, industry, and 
universities. The objectives of the workshop were to review the diverse activities in 
hierarchical approach to material modeling from continuum to atomistics; applications of 
multiscale modeling to advanced and improved material synthesis; defects, dislocations, 
and material deformation; fracture and friction; thin-film growth; characterization at nano 
and micro scales; and, verification and validation of numerical simulations, and to 
identify their potential for future aerospace systems. 
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INTRODUCTION 


In recent years intense effort has been devoted to the modeling and simulations of 
physical phenomena occurring on a vast range of spatial and temporal length scales. 
Multiscale modeling and simulations are being used in fields as diverse as materials 
science, nano/microelectronics, environmental remediation, safe stewardship of nuclear 
weapons and biotechnology. This endeavor has prompted the development of multiscale 
modeling and simulation strategies, interrogative visualization facilities, as well as 
versatile and powerful software systems. 

An attempt is made in this overview to give broad definitions to the terms and to set 
the stage for the succeeding presentations. The presentation is divided into four parts 
(Figure 1). In the first part, examples of future aerospace systems are given, along with 
some of their major characteristics and the enabling technologies. The second part 
provides a brief overview of some aspects of multiscale modeling and simulation as one 
of the enabling technologies for future aerospace systems. The third part describes the 
future research and learning environments required for the realization of the full potential 
of multiscale modeling and simulation. The fourth part lists the objectives of the 
workshop and some of the sources of information on multiscale modeling and simulation 


* Future Aerospace Systems 
and Enabling Technologies 

* Multiscale Modeling and 
Simulation 

* Future Research and 
Learning Environments 
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EXAMPLES OF FUTURE AEROSPACE SYSTEMS AND SOME OF THEIR 

CHARACTERISTICS 


The realization of NASA’s ambitious goals in aeronautics and space with the current 
national budget constraints will require new kinds of aerospace systems and missions that 
use novel technologies and manage risk in new ways. Future aerospace systems must be 
autonomous, evolvable, resilient, and highly distributed. Two examples are given in 
Figure 2. The first is a biologically inspired aircraft with self-healing wings that flex and 
react like living organisms. It is built of a multifunctional material with fully integrated 
sensing and actuation, and unprecedented levels of aerodynamic efficiencies and aircraft 
control. The second is an integrated human-robotic outpost, with biologically inspired 
robots. The robots could enhance the astronaut’s capabilities to do large-scale mapping, 
detailed exploration of regions of interest, and automated sampling of rocks and soil. 
They could enhance the safety of the astronauts by alerting them to mistakes before they 
are made, and letting them know when they are showing signs of fatigue, even if they are 
not aware of it. 
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ENABLING TECHNOLOGIES FOR FUTURE AEROSPACE SYSTEMS 


The characteristics of future aerospace systems identified in Figure 3 are highly 
coupled and require the synergistic coupling of the revolutionary and other leading-edge 
technologies listed in Figure 3. The three revolutionary technologies are nanotechnology, 
biotechnology and information/knowledge technology. The other leading-edge 
technologies are high-productivity computing; high-capacity communication; multiscale 
modeling, simulation and visualization; intelligent software agents; human performance, 
and human-computer symbiosis. 
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Figure 3 
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DEFINITION AND GOAL OF MULTISCALE MODELING 


The term multiscale modeling has been used in reference to various activities for 
exploiting insights arising either from distinct methodologies, or from attempts to 
incorporate multiple mechanisms, occurring at disparate spatial and/or temporal scales, 
into the same modeling paradigm. The overall goal of multiscale modeling is to predict 
the response of complex systems across all relevant spatial and temporal scales (Figure 

4 ). 
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SOME OF THE APPLICATIONS OF MULTISCALE 
MODELING AND SIMULATION 


5 ): 


Applications of multiscale modeling and simulation methodologies include (Figure 


1. Material synthesis and development 

2. Material characterization and degradation (e.g., damage, failure and plasticity) 

3. Phase transformations in materials (e.g., recrystallization) 

4. Turbulence modeling 

5. Atmospheric modeling 

6. Modeling of biological systems 

7. Systems on a chip 

8. Signal processing and analysis 



Figure 5 
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APPROACHES FOR MULTISCLAE MODELING 


Three categories of approaches can be identified for multiscale modeling, namely 
(Figure 6). 

■ Parameters passing, hierarchical modeling approaches in which a hierarchy of 
approaches and mathematical/computational models with different physical levels 
of description is pieced together. The output of the smaller-scale models is used 
as input for the larger-scale models. 

■ Handshaking approaches in which several computational approaches (and 
mathematical models) are used concurrently and use some sort of handshaking 
procedure to communicate 

■ Approaches that host more than on physical level of description in the same 
mathematical model. Example is provided by the Lattice-Boltzman equation 
which incorporates atomistic, kinetic, and continuum fluid description. 

Examples of the three categories are given subsequently. 
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EXAMPLES OF HIERARCHICAL MATERIAL MODELING 


Hierarchical material modeling has been used to predict and explain the mechanical 
properties of metals at dimensions ranging from a fraction of a nanometer to meters. The 
focus is on four major length scales - the atomic (nanoscopic) scale (nanometers), the 
microscale (micrometers), and mesoscale (millimeters), and the macroscale (centimeters 
and above). Fundamental physical and mathematical principles are rigorously applied to 
the modeling at each scale, and then data are then passed to the next scale up (Figure 7). 

For nanocomposite structures, six different levels can be identified in the hierarchy; 
namely, nanotubes and nanopolymers, fibers and matrix, ply, laminate, subcomponent 
and component. The techniques/theories used to go from each scale to the next are 
identified in Figure 7. In the hierarchical sensitivity analysis an opposite path is taken 
from that of the hierarchical modeling. The sensitivity coefficients with respect to each 
of the constituent and effective ply properties can be expressed as a linear combination of 
the sensitivity coefficients with respect to the laminate stiffnesses. 



Figure 7 
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EXAMPLES OF HANDSHAKING APPROACHES 


Two examples are given in Figure 8 for the multiscale modeling approaches that 
attempt to link several computational approaches in a combined model. The first is 
dynamic fracture analysis and the second is functionalized AFM tip (molecular robotics). 
In both examples, electronic structure model (quantum mechanics) is combined with a 
molecular dynamics model, which in turn is embedded into a continuum model 
(discretized) by finite elements. 



Molecular Robotics: Functionalized AFM Tip 



Figure 8 
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HIERARCHY OF SCALES IN FRACTURE 


The hierarchy of scale that may determine the response of a crack in a relatively 
ductile polycrystalline metal are shown in Figure 9. Zooming in to the near-tip region, 
the plastic zone governed by the macroscopic plastic flow is observed. Inside that zone 
the material is polycrystalline, with the plastic deformation being different in different 
grains. Closer to the tip of the crack plastic deformation occurs via localized slip on 
particular slip planes, caused by the motion of individual dislocations. The atomic 
arrangement on either side of the crack is shown in the bottom left part of the figure. 
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LENGHTH AND TIME SCALES IN HIERARCHICAL MODELING 


In computationally driven material development, a hierarchy of material models are 
used, associated with phenomena occurring at length scales from 10' 1IJ to lm and time 
scales from 10' 15 second to years (Figure 10). The disciplines involved in these scales 
include quantum mechanics; molecular dynamics; nanomechanics; micro, meso and 
macromechanics; process simulation and; engineering design. However, many gaps still 
exist in the hierarchy of models, and to date, no rational way exists to integrate these 
models and to couple them in order to relate the phenomena at the very small length 
scales with the macroscopic behavior. 



Figure 10 
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LATTICE-BOLTZMANN METHOD 


The Lattice-Boltzmann method (LBM) is an example of the modeling approaches that 
host more than one physical level of description (Figure 11). Unlike conventional 
methods which solve the discretized macroscopic Navier-Stokes equations, the LBM is 
based on microscopic particle models and mesoscopic kinetic equations. It is a derivative 
of the lattice gas automata methods, and is especially useful for modeling interfacial 
dynamics, flows over porous media, and multiphase flows. 



Length Scale 



Figure 1 1 
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QUANTUM MECHANICS 


Quantum mechanics is a means to understand and predict the interactions between 
atoms and molecules and to model chemical reactions at that scale (Figure 12). It uses 
models based on electronic structure. The solution of Schrodinger equation provides 
electronic wave functions. Other properties are obtained from these functions. 

The practical time and length scales are often a million times too small for industrial 
design. To overcome this difficulty, a hierarchy of methodologies is used with quantum 
mechanics at the foundation. 



Study of electronic structure of solids, surfaces or 


molecular scale 


* Predicts the interactions between atoms and molecules 


Hartree-Fock Approximation 


Figure 12 
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MOLECULAR DYNAMICS 


Molecular dynamics is a means to study matter at the atomic level and to predict the 
static and dynamic properties from the underlying interactions between the molecules 
(Figure 13). It allows the prediction of systems 100 to 10,000 larger than quantum 
mechanics for periods 1,000 to 1,000,000 times longer. 

To go from quantum mechanics to molecular dynamics requires averaging over the 
electrons to obtain spring constants, discrete charges and van der Waals parameters. 

Recently, a team of researchers from Lawerence Livermore National Laboratory, 
IBM, and the Max Planek Institute for Metal Research in Germany were able to simulate 
the interactions of one billion atoms and discover cracks traveling faster than the speed of 
sound (767 mph in air at sea level and room temperature). They used half the 8,192 
processors of the ASCI white computer for five full days (about million processor hours). 



Figure 13 
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NANOMECHANICS 


Nanomechanics deals with mechanics problems associated with modeling, design, 
fabrication and application of 3-D structures and systems with submicron (nm-scale) 
dimensions (Figure 14). Nanoscale systems have a number of interesting features which 
distinguish them from micro- and marco-scale systems. 



Figure 14 


16 




CONTINUUM MECHANICS VS. NANOMECHANICS 


Continuum mechanics is used to predict the phenomena described by uniform 
collective behavior of atoms (for example, elastic response of materials). By contrast, 
nanomechanics is used to predict the phenomena described by dramatic changes in the 
state of few atoms (for example, corrosion and fatigue). To accomplish this, the atoms 
whose state is affected need to be identified, their behavior at the nm length scale 
modeled and connected to the solid’s macroscopic response (Figure 15). 



Figure 15 
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KEY COMPONENTS OF ADVANCED 
MULTSCALE RESEARCH ENVIRONMENT 


Multiscale modeling is a unifying paradigm to enable integration of science and 
engineering. It allows rigorous correlation of different science and engineering models, 
representations, languages and metrics. It overcomes engineering process fragmentation. 

The realization of the strategic value of multiscale modeling requires an advanced 
research environment that links diverse teams of scientists, engineers and technologists. 
The essential components of the environment can be grouped into three categories: 
intelligent tools and facilities, nontraditional methods, and advanced interfaces (Figure 
16). The three categories are described subsequently. 



Figure 16 
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INTELLIGENT TOOLS AND LACILITIES 


These include high fidelity - rapid modeling, simulation and interrogative 
visualization tools, synthetic immersive environments; close coupling of simulations and 
experiments; computer simulation of physical experiments and remote control of these 
experiments. In all of these tools, extensive use should be made of intelligent software 
agents and information technology (Figure 17). 
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Figure 17 
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NONTRADITIONAL METHODS 


These include novel mathematical tools for multiscale data representation (e.g., 
wavelets, ridgelets and curvelets) and multiresolution methods; verification and 
validation strategies; and nondeterministic approaches for handling various types of 
uncertainties in geometry, material properties, boundary conditions, loading and 
operational environment (Figure 18). 



Figure 18 


20 





MULTISCALE DATA PRESENTATION 


In recent years a number of novel concepts related to multiscale data representation 
and wavelet bases have appeared. Wavelets are functions of compact support - small 
waves. They are well-suited for describing local phenomena (e.g., point of zero- 
dimensional singularities). They oscillate and have zero net area. 

In dimensions greater than one ridgelets and curveiets are used for representing sharp 
variations, or singularities. Ridgelets are used for linelike and planelike phenomena and 
curveiets are used for representing sharp variations across curves (Figure 19). 
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MULTIRESOLUTION METHODS 


In typical multiresolution methods a sequence of system descriptions at a sequence of 
representative scales are constructed, combining local processing at each scale with 
various inter-scale interactions. For example, the evolving solution at each scale is 
typically used to construct the equations on coarser scales and to accelerate the solution 
on finer scales, hi this way, large-scale changes are effectively performed on coarse 
grids, based on information previously gathered from finer grids (Figure 20). 




Figure 20 
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VERIFICATION AND VALIDATION OF NUMERICAL SIMULATIONS 


Quantifying the level of confidence, or reliability and accuracy of numerical 
simulations has recently received increased levels of attention in research and 
engineering applications. During the past few years, new technology development 
concepts and terminology have arisen. Terminology such as virtual prototyping and 
virtual testing is now being used to describe computer simulation for design, evaluation 
and testing of new engineering systems. 

The two major phases of modeling and simulation of an engineering system are 
depicted in Figure 21. The first phase involves developing a conceptual and 
mathematical model of the system. The second phase involves discretization of the 
mathematical model, computer implementation, numerical solution and representation or 
visualization of the solution. In each of these phases there are uncertainties, variabilities 
and errors. 

Verification and validation are the primary methods for building and quantifying 
confidence in numerical simulations. Verification is the process of determining that a 
model implementation accurately represents the conceptual/mathematical model and the 
solution to the model. Correct answer is provided by highly accurate solutions. 
Validation is the process of determining the degree to which a model is an accurate 
representation of the real system. Correct answer is provided by experimental data. 



Figure 21 
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TYPES OF UNCERTAINITES 


A number of different uncertainty representations and classification have been 
proposed. Among these classifications are the three-type classification - statistical, 
model, and fundamental uncertainties; the two type classification - uncertainty of 
information and uncertainty of the reasoning process, and the six-type classification (see 
Figure 22): 

■ Probabilistic uncertainty, which arises due to chance or randomness, 

■ Fuzzy uncertainty due to linguistic imprecision (e.g., set boundaries are not 
sharply defined) 

■ Model uncertainty which is attributed to lack of information about the model 
characteristics, 

■ Uncertainty due to limited (fragmentary) information available about the system 
(e.g., in the early stage of the design process), 

■ Resolutional uncertainty which is attributed to limitation of resolution (e.g., 
sensor resolution), and 

■ Ambiguity (e.g., one to many relations) 

Accounting for uncertainties in the design of engineering systems involves 
representation and propagation of uncertainty, from model parameters to estimating the 
performance; bounding uncertainties; model updating and validation; and assessment of 
the impact of uncertainties on the reliability and certification of the system. 




Figure 22 
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PRINCIPLE OF COMPLEXITY 


One of the important consequences of uncertainty is its effect on precision. Three 
types of models can be identified depending on the complexity and the precision, namely: 
mathematical models, model-free methods, and fuzzy systems. In a typical complex 
system a combination of the three should be used. As the uncertainty and/or complexity 
of an engineering system increases, the ability to predict its response diminishes, until a 
threshold is reached beyond which precision and relevance become almost mutually 
exclusive. Consider, for example, numerical simulations in which sophisticated 
computational models are used for predicting the response, performance, and reliability 
of the engineering system, but the system parameters are little more than guesses. Such 
simulations can be characterized as Correct but Irrelevant Computation (CBIC); that is, 
forcing precision where it is not possible (Figure 23). 



Figure 23 
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NONDETERMINISTIC ANALYSIS APPROACHES 


Three general approaches can be used for the analysis of systems with uncertainties; 
namely (Figure 24): probabilistic methods for random processes; fuzzy sets; and set 
theoretical or antioptimization methods. The domain of application of each of these 
techniques is identified in Figure 24. 
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ADVANCED HUMAN/COMPUTER INTERFACES 


Although the WIMP (windows, icons, menus, and pointing devices) paradigm has 
provided a stable global interface, it will not scale to match the myriad form factors and 
uses of platforms in the future collaborative distributed environment. Perceptual user 
interfaces (PUI’s) are likely to meet those needs. PUI’s integrate perceptive, multimodal 
and multimedia interfaces to bring human capabilities to bear on creating more natural 
intuitive interfaces. They enable multiple styles of interactions, such as speech only, 
speech and gesture, vision, and synthetic sound, each of which may be appropriate in 
different applications (Figure 25). These new technologies will enable broad uses of 
computers as assistants or agents that will interact in more human-like ways. 



* Integrates perceptive, multimodal and multimedia 
interfaces to bring human capabilities to bear on 
creating more natural and intuitive interfaces 

* Enables multiple styles of interactions and broad 
uses of computers as assistants 


Figure 25 


27 






ADVANCED LEARNING ENVIRONMENTS 


In order to meet the life-long learning demands of the future and broaden the 
awareness among the researchers and engineers of nondeterministic approaches, three 
categories of learning environments are needed; namely, expert-led group learning 
environment; self-paced individual learning environment; and collaborative learning 
environment (Figure 26). The three environments, in combination, can reduce the time 
and cost of learning, as well as sustain and increase worker competencies in high tech 
organizations. 

The human instructors in these environments will serve many roles, including 
inspiring, motivating, observing, evaluating, and steering the learners, both individually 
and in distributed teams. 



Figure 26 
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EXPERT-LED LEARNING ENVIRONMENT 


The human instructors in expert-led distributed learning in a virtual environment 
serve as coaches, guides, facilitators, and course managers. Their presentations focus on 
a broad overview of the topic and its diverse applications (Figure 27), and end with more 
penetrating, what-if questions that can enhance the critical thinking and creativity of the 
learners. Elaborate visualization and multimedia facilities are used in the presentations. 
Routine instructional and training tasks are relegated to the self-paced individual learning 
environment. 
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SELF-PACED LEARNING ENVIRONMENT 


The individual learning environment engages the learner and provides a high degree 
of tailored interactivity. It can be used for self-paced instruction of routine material not 
covered in the lecture. Using virtual instructors assigned by the human instructors can 
enhance such instruction. It can be used to study the physical phenomena occurring at 
different length scales using advanced visualization, multimedia and multisensory 
immersive facilities. The individual learning environment can serve to carry out 
numerical and virtual experiments - computer simulation of physical experiments (Figure 
28). 



Figure 28 
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SELF-PACED LEARNING ENVIRONMENT 
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COLLABORATIVE / DISTRIBUTED LEARNING ENVIRONMENT 


Collaborative learning environments teach teamwork and group problem solving. 
Instructors and learners can be geographically dispersed. Eventually, they can be brought 
together through immersive telepresence facilities to share their experiences in highly 
heterogeneous environments involving different computing platforms, software, and 
other facilities, and they will be able to work together to design complex engineering 
systems beyond what is traditionally done in academic settings. Because participants can 
be virtually collocated without leaving their industry and government laboratories, 
collaborative learning environments can enable the formation of learning networks 
linking universities, industry and government labs. The ultimate goal of these learning 
facilities is to create an intellectual environment where academic and experiential 
learning are effectively and efficiently co-mingled. In such an environment, academic 
rigor is learned in concert with professional job performance, and academic complexities 
are addressed within the industrial concern (Figure 29). 



Figure 29 
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COLLABORATIVE / DISTRIBUTED LEARNING ENVIRONMENT 
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VIRTUAL CLASSROOM 


Online training and virtual classrooms are typically used to provide learning 
environments with custom self-instruction, flexible tutorial support, and choice of both 
the place and time of learning. Three categories of facilities are used in these 
environments; namely: instruction, including multimedia lectures, links to other 
resources and tools for searching, browsing, and using archived knowledge; 
communication, including email, UseNet, chat centers, video and Internet conferencing; 
and course management and performance evaluation (Figure 30). 
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MULTISCALE MODELING AND SIMLUATION NETWORK 


The realization of the full potential of multiscale modeling and simulation in the 
design, and development of future complex systems requires, among other things, the 
establishment of research and learning networks. The networks connect diverse, 
geographically dispersed teams from NASA, other government labs, university consortia, 
industry, technology providers, and professional societies (Figure 31). 



Figure 31 


35 




EVOLUTION OF NEW TECHNOLOGY 


The multiscale modeling, simulation and visualization, and their associated 
technologies, as any other technology, have gone through three phases. The first is that 
of naive euphoria - unrealistic expectations resulting from overreaction to immature 
technology. The second is cynicism , or frustration associated with unmet expectation. 
The third is that of realistic expectations - gradually realizing the true benefits from the 
technologies (Figure 32). 



Figure 32 
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OBJECTIVE AND FORMAT OF WORKSHOP 


The objectives of the workshop are to: a) provide a broad overview of the diverse 
activities related to multiscale modeling, simulation and visualization, and b) identify the 
potential of these technologies to future aerospace systems (Figure 33). The nineteen 
presentations illuminate some of the key issues in multiscale modeling and simulation, 
particularly in the area of material modeling and synthesis, and provide fresh ideas for 
future research and development. 


Objectives 

* Overview of diverse activities related to multiscale modeling, 
simulation and visualization 

* Identity potential for future aerospace systems 

Format 

* 1 9 presentations, 7 sessions 

* Brainstorming session 


* Printed (NASA CP) 


Figure 33 
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ISSUES RELATED TO MULTISCALE MODELING, 
SIMULATION AND VISUALIZATION 


The issues related to, and the pacing items for, multiscale modeling, simulation and 
visualization can be grouped into five major categories, and are identified in Figure 34. 


General 

* Rote and potential of multiscale modeling and simulation in 

engineering design - current barriers and challenges 

* Synergistic coupling of multiscale simulations and physical 

experiments - mathematical models and physical experiments 
evolving together, with each needing the other 

* Role of virtual and hybrid physical/virtual tests 

Multiscale Linking /Modeling 

* Accurate mathematical models describing heterogeneous systems, 

multiphysics/muftiphenomena of interest at disparate spatial/ 
temporal scales 

» Account for the interconnectivity of the essentia! phenomena 
occurring at the spectrum of the length scales (including, detailed 
molecular level interactions, and preserving the macroscopic 
conservation principles). 

« Model robustness - emergent behaviors insensitive to model details 
and starting conditions 

Physical Experiments/Design of Experiments 

« Experimental techniques and facilities needed to identify important 
phenomena, particularly at the small length and time scales 
« Experimental data needed to initialize and/or validate simulations 

* Design of validation experiments (measuring reality) 

* Definitions of feature fidelity metrics 

* Variability: Replication of experiments 

* Estimating and accounting for experimental uncertainties 

Multiscale Simulations 

« Seamless integration of different length and time scales, and creation 
of unified multiscaie simulation tools 
« Intrinsic limits to simulation (computability) 

* Scaling algorithms to increased model size and resolution 

* Building intelligence into the simulation 

* Benchmark simulations at various length scales for establishing 

reliability of simulations 

Visualization 

* Role of interrogative visualization (e.g., feature extraction), 
computational steering and inverse steering 


Figure 34 
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INFORMATION ON MULTISCALE MODELING, SIMULATION 
AND VISUALIZATION 


Extensive literature now exists on multiscale modeling and simulation. A short list of 
books, monographs, conference proceedings, survey papers and websites is given 
subsequently. 

Books. Monographs and Conference Proceedings 


■ Phillips, R., “ Crystals , Defects and Micro structures: Modeling Across Scales’ ’, 
Cambridge University Press, Cambridge, UK. 2001. pp. 780. 

■ Kubin, L. P., Selinger, R. L., Bassani, J. L., and Cho, K. (editors), “ Multiscale 
Modeling of Materials - 2000 ”, Proceeding of a symposium held Nov. 27 - Dec. 
1, 2000, Boston, MA, Materials Research Society, Volume 653, 2001. 

■ Murayama, Y., “ Mesoscopic Systems - Fundamentals and Applications'”, John 
Wiley & Sons, Inc. Chichester, United Kingdom, 2001. 

■ Ohno, K., Isfarjani, K., and Kawazoe, Y., “Computational Materials Science - 
from Ab Initio to Monte Carlo Methods'”, Springer- Verlag, Berlin, Germany, 

1999. 

■ Cronin, J. and O’Malley, R.E. (editors), “Analyzing Multi scale Phenomena Using 
Singular Perturbation Methods”, Proceedings of Symposia in Applied 
Mathematics, Volume 56, American Mathematical Society, 1999. 

■ Bulatov, Vasily, V., Diaz de la Rubia, T., Phillips, R., Kaxiras, E., and Ghomem, 
N. (editors), “Multiscale Modelling of Materials” , Proceedings of a symposium 
held Nov. 30 - Dec. 3, 1998, Boston, MA, Materials Research Society, Volume 
538, 1998. 

■ Raabe, D., “ Computational Materials Science - The Simulation of Materials 
Microstructures and Properties”, Wiley-VCH Verlag GmBH. D-69469 
Weinheim (Federal Republic of Germany), 1998. 

■ Golden, K.M., Grimmet, G.R., James, R.D., and Milton, G.W. (editors), 

“ Mathematics of Multiscale Materials” , Springer- Verlag New York, Inc., 1998. 

■ Dahmen, W., Kurdilla, A.J., and Oswald, P. (editors), “Multiscale Wavelet 
Methods for Partial Differential Equations (Wavelet Analysis and Its 
Applications, V. 6), Academic Press, 1997. 

■ Hiller, J. R., Johnston, I. D., and Styer, D. F., “Quantum Mechanics Simulations - 
The Consortium for Upper-Level Physics Software ”, John Wiley & Sons, Inc., 
New York, 1995. 

■ Rapaport, D., “ The Art of Molecular Dynamics Simulation” , Cambridge 
University Press, Cambridge, United Kingdom, 1995. 

■ Cohen, A., “ Wavelets and Multi scale Signal Processing” , CRC Press, 1995. 
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Journal Articles 


■ Needleman, A., and Van der Giessen, E., “Micromechanics of Fracture: 
Connecting Physics to Engineering” , MRS Bulletin, Materials Research Society, 
Volume 26, No 3, March 2001. 

■ Diaz de la Rubia, T. D., and Bulatov, V. V., (guest editors) “Materials Research 
by Means of Multiscale Computer Simulation”, MRS Bulletin, Materials 
Research Society, pp. 169 - 170, March 2001. 

■ Tolies, W. M., “ Self-Assembled Materials ”, Materials Challenges for the Next 
Century, MRS Bulletin, Materials Research Society, October 2000, pp. 36 - 38. 

■ Hafner, J., “ Atomic-Scale Computational Materials Science ”, Acta Materiala, 
Volume 48 (2000), pp. 71 - 92. 

■ Needleman, A., “Computational Mechanics at the Mesoscale” , Acta Materiala, 
Volume 48 (2000), pp. 105 - 124. 

■ Wild, A., “ Multi-discipline , Multi-scale modeling of Microsystems: An 
Overview ”, Proceeding of the IEEE International Caracas Conference on Devices, 
Circuits and Systems, ICCDCS, 2000, Piscataway, NJ, USA, pp. Cl 10-1 - Cl 10- 
5. 

■ Lemarchand, C., Chaboche, J. L., Devincre, B., and Kubin, L. P., “ Multiscale 
Modeling of Plastic Deformation”, Journal de Physique IV, Proceedings of the 
1998 3 ,d European Mechanics of Materials Conference on Mechanics and Multi — 
physics Processes in Solids: Experiments, Modeling, Applications (EUROMECH 
- MECAMAT ’98), Oxford, United Kingdom, 1 999. 

Websites 


■ “ Predicting Material Behavior from the Atomic Level Up ”, Multiscale Modeling ” 

http : //vvvvw. llnl . gov/str/Mori art v . him] . 

■ “ Exploring the Fundamental Limits of Simulation” 

http : //www. llnl, gov/ s tr/MarchO i /Gal li . h tin 

■ “Multiscale Materials Modeling Presentations” 

http : L L www. tc . com ell . edu/Research/Multi seal e/Presentations/ 

■ “Visual Quantum Mechanics” 

http://phvs. edue. ksu. edu/ 
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CONCURRENT COUPLING OF LENGTH SCALES 


This work on concurrent coupling of length scales was started by Farid Abraham at 
IBM Almaden and Jeremy Broughton (then at NRL), together with Prof. Tim Kaxiras at 
Harvard and myself. The work has been supported by the Office of Naval Research and 
the Naval Research Laboratory, and computer time was provided by the DOD HPCMO 
Challenge Program. 


CONCURRENT COUPLING OF LENGTH 

SCALES 


Noam Bernstein 

Center for Computational Materials Science 
Nava! Research Laboratory Washington, DC 


Collaborators: 

F. F. Abraham 
E. Kaxiras 
D. W. Hess 

Support: 

NRL and ONR 
DOD HPCMO 
Challenge Program 



Figure 1 
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OVERVIEW 


The talk begins with an overview of the coupling of length scales (CLS) concept. In 
particular, I discuss the differences between “sequential” and “concurrent” coupling of 
length scales approaches. After a brief introduction to the relevant issues in fracture, I 
describe our first approach to the CLS method. I describe the method we developed, and 
some results and lessons we learned from that work. I then describe our current CLS 
approach for simulating fracture. After presenting the method and results, I present some 
future directions and conclusions. 


OVERVIEW 

Coupling of length scales 
sequential vs. concurrent 

Fracture: the canonical multiscale materials problem 

CLS first approach 
method 
results 

lessons learned 

CLS: new approach 
method 
results 

Directions for future 
Conclusions 


Figure 2 
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COUPLING OF LENGTH SCALES 


Many material properties are controlled by processes that occur on a wide range of 
length scales. Each of these processes is most appropriately (accurately and efficiently) 
described by different physical models. 

One prototypical example is fracture. On the shortest length scales there are a few 
hundred atoms near the crack tip where breaking of interatomic bonds is significant. For 
this region a classical description of the nuclei is sufficient, but the electrons should be 
described quantum-mechanically for the greatest accuracy. Farther from the crack tip 
there are tens of thousands of atoms where bonds are highly strained but not broken. In 
this region empirical potentials (EPs) do a good job of describing the interactions 
between the atoms. The elastic deformation that drives fracture can extend over much 
larger, macroscopic length scales, and in the vast majority of this region is very 
accurately and efficiently described by continuum mechanics. 

How to best describe a system with all of these different processes and regions is still 
an open question. 


COUPLING OF LENGTH SCALES 


Material properties controlled by processes over many length scales 


Different processes best described by different models 


Example: fracture 

• Short scale - io 2 atoms 

- breaking atomic bonds 

- classical nuclei, quantum-mechanical bonds 

• Medium scale - 10 4 atoms 

- highly strained bonds 

- classical nuclei, empirical interactions 

• Long scale - IQ 23 atoms 

- elastic deformation 

- continuum mechanics 

How can you treat all of these aspects? 



Figure 3 
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SEQUENTIAL COUPLING 


Sequential coupling is an approach where an accurate method is used to simulate a 
number of configurations in a small system. From these simulations some parameters for 
a faster but more empirical method are extracted. One popular example is kinetic monte 
carlo for surface growth, where barriers to adatom motion can be computed with density- 
functional theory. Another is discrete dislocation dynamics, where the mobilities are 
computed using atomistic simulations. 

In any implementation some general features are needed for the sequential approach 
to apply. One is a separation of time scales. If the essential events of the faster method 
occur on the same time scales as the details of the processes computed using the slower 
system, the two can’t be separated. The other major requirement is that there exists a 
faster and more empirical theory to use, and that it has a reasonable number of parameters 
than can be practically computed using the more accurate method. 
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CONCURRENT COUPLING 


An alternative to the sequential approach is concurrent coupling, where a large 
system is simulated using a fast method with limited applicability, and where needed a 
slower, more accurate method is used. Two good examples where this approach is 
important are fracture and friction. In fracture the crack tip bond breaking can be 
described with a quantum-mechanical model of bonding, while the rest of the sample is 
described with empirical potentials. In friction it might be necessary to describe surface- 
surface interaction using quantum-chemical approaches while using continuum elasticity 
to simulate the contact forces. 

The essential requirement for the concurrent approach is that the slower method is 
only required to describe a localized region. The main importance of the concurrent 
approach is that it can be applied to dynamics, and to systems where the boundary 
conditions on the region that must be described accurately are changing with time. In 
these situations information is flowing in both directions between the regions described 
by the different methods. This makes it difficult to summarize the effect of one on the 
other in terms of a boundary conditions or a few parameters. 


CONCURRENT COUPLING 

Concurrent: large system with a fast method, limited applicability 
where needed, use a slower, more accurate method 

• Fracture: QM for crock tip. bulk sample with empirical potentials 

• Friction: QC for surface Interaction, elasticity for contact forces 
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Mechanics 
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* localized region for slower method 


Especially Important for dynamics, changing boundary conditions 


Figure 5 
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FRACTURE 


Fracture is a phenomenon of great technological importance. It would be very useful 
to be able to predict fracture properties, ideal from first-principles. It would also be useful 
to be able to understand fracture properties well enough to control them. Fracture is also 
interesting for purely scientific reasons. Why is it that materials break in the way they do, 
some through brittle fracture, others through plasticity and ductile fracture? 

Silicon in particular is an interesting example of the transition between brittle and 
ductile fracture. Over a few degrees temperature variation its stress to fracture changes by 
almost a factor of two. 1 For this reason it has become a model system for the brittle to 
ductile transition. 

1 P. B. Hirsch and S. G. Roberts, “The brittle-ductile transition in silicon,” Phil. Mag. 
A, Vol. 64, 1991, pp. 55-80. 
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ORIGINAL CLS METHOD 


The original CLS method 2 ’ 3 used a well defined Hamiltonian to express the total 
energy of the system. This energy consisted of contributions from three regions: 
continuum elasticity, empirical potentials, and tight-binding (TB). The first was 
evaluated using 2-D finite elements. The second used the Stillinger- Weber (SW) 
interatomic potential for Si. 4 The third used the Si TB model of Bernstein and Kaxiras 5 
with an approximate solution for the forces on the atoms. The dynamics for each region 
were advanced forward in time in sync. 

2 F. F. Abraham, J. Q. Broughton, N. Bernstein, and E. Kaxiras, “Spanning the 
continuum to quantum length scales in a dynamic simulation of brittle fracture,” 
Europhys. Lett., Vol. 44, 1998, pp. 783-787. 

3 J. Q. Broughton, F. F. Abraham, N. Bernstein, and E. Kaxiras, “Concurrent coupling 
of length scales: Methodology and application,” Phys. Rev. B, Vol. 60, 1999, pp. 2391- 
2403. 

4 F. H. Stillinger and T. A. Weber, “Computer simulation of local order in condensed 
phases of silicon,” Phys. Rev. B, Vol. 31, 1985, pp. 5262-5271. 

5 N. Bernstein and E. Kaxiras, “Nonorthogonal tight-binding Hamiltonians for defects 
and interfaces in silicon,” Phys. Rev. B, Vol. 56, 1997, pp. 10488-10496. 


ORIGINAL CLS METHOD 



Conservative H a i ni i ton ia n : 

///,* = E «fc+Wjjp+ E 11 m 

FE regions TB regions 

• Vf fj Continuum elasticity (Hooke's law) 

• \) n Stltlinger-Weber interatomic potential 

• • j-j> Tight-binding (Bernstein and Kaxiras) 

Dynamics for each region in sync 



Figure 7 


49 





ORIGINAL CLS RESULTS 


The picture shows a visualization of the crack tip after substantial propagation using 
the original CLS method without the tight-binding region. The whole sample is an 1 1 A 
thick, periodically repeated system about 4000 A by 3600 A on a side. There are 1.5 x 
10 6 empirical potential atoms, and another three 3 x 10 5 finite element nodes representing 
about 7.7 x 10 6 atoms. 

In the image, color coded by the local energy of each atom, there are elastic waves 
emanating from the crack, as well as voids and amorphous tendrils near the crack. The 
amorphous tendrils are left behind rapidly propagating dislocations. 

This type of morphology is indicative of fracture that is not at all brittle. The energy 
release rate, a measure of the elastic energy that’s driving the fracture process is about 
130 J/ 3 , many times larger than the experimental results of about 2.5 J/nT. 7 

7 J. A. Hauch, D. Holland, M. P. Marder, and H. L. Swinney, “Dynamic fracture in 
single crystal silicon,” Phys. Rev. Lett., Vol. 82, 1999, pp. 3823-3826. 


ORIGINAL CLS RESULTS 



4000 A x 3600 A x 11 A 
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Elastic waves 
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LESSONS FROM ORIGINAL CLS 


The original CLS approach combined continuum mechanics, empirical potentials, and 
tight-binding. The results showed that the simulations were stable, and efficient enough 
to run for long times. However, the empirical potentials were not accurate enough to 
describe the bond breaking process, nor were the approximations used to evaluate the TB 
forces. There was also a wide range of relevant time scales between the finite element 
and tight-binding regions. One conclusion is that it may not be obvious apriori what will 
prove to be a sufficiently accurate method. The other is that the simulation duration is 
limited by the computational demands of the TB region. 

A back-of-the-envelope estimate of the time scales is simple: The length of the run is 
limited by the most accurate method, TB, which take about 10 s for a 1 fs time step. One 
day of runtime is about 1000 time steps, or about 1 ps. There is a characteristic velocity, 
the speed of sound, which is of order 4 km/s, giving a length scale of about 40 A. A 
region on this length scale can almost always be simulated efficiently enough 
atomistically, and the complications involved in coupling atomistics to continuum 
mechanics are not needed. 

The inaccuracies in the approximations used for the TB, together with the separation 
of time scales, led us to modify the approach without finite elements and a more accurate 
solution of the TB forces. 


LESSONS FROM ORIGINAL CLS 

continuum mechanics, empirical potentials, tight, binding 
Results: 

♦ Stable, long running simulations 

♦ Empirical potential, tight-binding solution not accurate enough 

♦ Wide range of time sales: FEM vs. TB 

Conclusions: 

♦ Sufficient accuracy may not he obvious a priori 

♦ Simulation duration limited by TB region 

Rut: is limited by most accurate method: TB 10 s/i n rime step 
t day of simulation =*• 1000 time stops ~ >s 

No FEM, improved TB 
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NEW CLS 


The new concurrent CLS approach used only atomistic methods to calculate the 
forces on atoms, which were then used to do molecular dynamics (MD) simulations. 
Most of the atoms in the fracture system were described using the EDIP empirical 
potential. 8 

The TB region was described with the same model as before 5 , but the forces were 
evaluated using a different method. 9 This method computed the single electron density 
matrix in real space by summing the Green’s function at a number of complex energies. 
The complex energies are computed from the poles of an approximate Fermi function. 
One advantage to this method is that it is possible to compute the density matrix in a 
localized portion of a large system. Another is that using the sparsity of the TB model 
Hamiltonian and overlap matrices, one can compute a sparse Green’s function matrix. 

8 J. F. Justo, M. Z. Bazant, E. Kaxiras, V. V. Bulatov, and S. Yip, “Interatomic 
potential for silicon defects and disordered phases,” Phys. Rev. B, Vol. 58, 1998, pp. 
2539-2550. 

9 N. Bernstein, “Linear scaling nonorthogonal tight-binding molecular dynamics for 
nonperiodic systems,” N. Bernstein, Europhys. Lett., Vol. 55, 2001, pp. 52-58. 


NEW CLS 

molecular dynamics - constant T, \ 

Empirical potentials: 

• EDIP (or Stllilnger-Weber) 

Tight, binding: 

• nonorthogonal, *p 3 . Bernstein and Kaxiras (or NRL-TB'} 

• density matrix from sum of Green’s functions (/(*,) 

• approximate Fermi function - small number of poles c. 

• linear scaling: 

compute sparse C(z) by numerically inverting -s' - //. 
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COUPLING AND TB 


An essential ingredient in the coupling scheme is computing the Green’s function, 
and therefore the density matrix, in a finite part of a larger system. To do this we compute 
the Green’s function using iterative inversion, applying a constraint to the Green’s 
function matrix elements at the boundary of the TB region. 

The mechanical coupling between the two regions is achieved simply. The EP atoms 
are included in the EP calculation, and are subject to forces from the EP calculation. The 
TB atoms are included in the TB calculation, and are subject to forces from the TB 
calculation. The boundary atoms are included in both calculations, but are subject only to 
EP forces. 


COUPLING EP AND TB 

G(z) In finite part of system 

Constrain C(z) matrix elements at boundary during matrix inversion 
Mechanical Coupling: 

* : included in EP calculation, forces from E P 

* TB atoms: included in TB calculation, forces from TB 

* boundary atoms: included in both calculations, force s from EP 



Figure 1 1 
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FRACTURE WITH TB 


Using this approach, we simulated fracture in silicon coupling the EDIP empirical 
potential to the TB model of Bernstein and Kaxiras. A seed crack was created in the 
system, and the sample was loaded in tension (mode I) by moving the atoms along the 
displacement field for the continuum elasticity solution for a finite crack in an infinite 
plate. 10 The entire system consists of about 5 x 10 4 atoms, with about 1000 atoms in the 
TB region. This results in a sample size of about 400 A by 250 A, in a 12 A thick, 
periodically repeated slab. The image below shows a 80 A by 65 A region around the 
crack tip. 

10 K. B. Broberg, Cracks and Fracture, Academic Press, San Diego, 1999. 


FRACTURE WITH TB 

Couple empirical potential (EDIP) and TB 
Continuum solution for fixed strain, top/bottom boundary fixed 
~50G00 EP atoms (1 proc.) ~1QG0 TB atoms (50 procs. ) 
400 A x 250 A x 12 A, (SO A x 65 A shown) 



Red: EP 
Blue: Boundary 
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BRITTLE FRACTURE 


The image on the left shows a visualization of the crack at the end of the MD 
simulation. The Si sample has failed through brittle fracture. The crack has propagated 
via interplanar cleavage, remaining atomically sharp and leaving behind a smooth crack 
face. 

On the right is a plot of the crack speed as a function of loading, measure through the 
energy release rate G. The vertical line is the Griffith criterion, the point at which the 
elastic energy reduction due to crack propagation is exactly equal to the surface energy of 
the increased exposed surface. The blue line indicates the results of the present 
simulation, in very good agreement with the black dots representing the experimental 
results of Hauch et al 1 The red line approximates the results for the empirical potentials, 
which fracture only at much higher loadings than experiment or the new CLS approach. 

The crack speed, which is about half of the Rayleigh wave speed, is similar in all sets 
of results. 


BRITTLE FRACTURE 
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DISCUSSION 


The results of the new CLS with an accurate calculation of the TB forces are clearly 
different from the EP results. However, it is unclear how the quantum-mechanical 
description of bonding makes a difference - what are the relevant fundamental material 
parameters that are different in the two descriptions? 


DISCUSSION 


EDIP EDIP+TB 



What’s different? 

I.e. what are the fundamental materials parameters that control the 
nature of fracture? 


Figure 14 
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ENERGIES 


One important class of approaches for predicting the fracture properties of a material 
is based on energetics. For example, the Rice criterion makes a quantitative prediction 
about whether a material will be brittle or ductile based on the ratio between the surface 
energy and the unstable stacking fault energy, which is the energy barrier to nucleating or 
propagating a dislocation. 1 1 

However, a comparison of this ratio for density-functional theory, 12 the TB model 
used here, and the two empirical potentials, 8 , does not show a systematic difference 
between the TB model and both of the EP models. Apparently the Rice criterion does not 
provide a good explanation for the differences between the EP simulation and the CLS 
method using the TB model. 

11 J. R. Rice, “Dislocation nucleation from a crack tip: an analysis based on the 
Peierls concept”, J. Mech. Phys. Solids, Vol. 40, 1992, pp. 239-271. 

12 E. Kaxiras and M. S. Duesbery, “Free energies of generalized stacking faults in Si 
and implications for the brittle-ductile transition,” Phys. Rev. Lett., Vol. 70, 1993, pp. 
3752-3755. 


ENERGIES . . . 

Energetics (Rice criterion): Compare 

• 7 s surface energy : make new crack surface 

* 7 us unstable stacking fault energy : make dislocations 
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. . . VS. FORCES 


An alternative to energy criteria is a criterion based on forces, independently 
proposed in various forms by Abraham and Marder. 

Forces depend on the way energy varies over a displacement distance. For 
dislocations, the distance is set by the lattice constant, which is the same for all 
reasonable models of a given material. For surfaces, on the other hand, the distance is set 
by two factors. In reality it is set by the physics of bonding. Any model of bonding makes 
some approximations that can greatly change this range. For example, first-principles 
methods in general involve much longer range interactions than empirical potentials, 
especially for covalently bonded materials such as Si. 

Since the surface energies for all the models are reasonable, the range of interaction 
ranges in the different models can greatly vary the forces that atoms feel as two surfaces 
are pulled apart. 


. . . VS. FORCES 


(Abraham, Marder) 

Force depends on energy /distance 


Dislocations (->•„„): distance set by lattice, same for ali models 
Surfaces (y s ): distance set by range of interactions in model 

♦ Determined by physics (covalent vs. Coulomb) 

• Restricted by model (DFT vs. empirical potentials) 



& 8 
& & & 
& & % 
^ & «c » 

f 9 

* ^ $ 


& # « II 

# # «r T 
» & & 

# & « 

$ « « 

« « # 


\ - 

# & & 

& 8 Sr 

# » ® 


# & 

« <£ 
# # & 
# « 


« & & % * 


A 

/ \ 


\ / 


\ / 
\ / 



Figure 16 


58 





RANGE: ENERGIES VS. FORCES 


The plot on the left is the energy required to rigidly pull apart two slabs of material as 
a function of the distance between them. This energy is normalized by the relevant 
comparable dislocation related energy, the unstable stacking fault energy. The plot on the 
right is the derivative of the energy as a function of distance, normalized by the Peierls 
stress which is the maximum derivative of the unstable stacking fault energy path. 

While the difference in ranges between the first-principles LDA, TB, and EP 
calculations is apparent, it does not show any systematic variations. The maximum force 
for the TB model is comparable to that of EDIP, and in fact much lower than for SW. If 
this force criterion were applicable, the TB simulation would look at least qualitatively 
like EDIP and unlike this is not the case, this type of force based criterion does not 
explain the differences either. 


RANGE: ENERGIES VS. FORCES 


Opening energy Opening stress 

(Normalized by %,*) (normalized by Peierls stress) 



vacuum thickness (A) vacuum thickness (A) 


Ranges differ, but not systematically 

Does not explain differences between EP and TB 

(LDA Peierls stress from Joos et 3i.) 


Figure 17 
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WHAT IS CAUSING THE DIFFERENCE? 


Energetics does not seem to explain the differences between the CLS and EP 
simulations, since the Rice criterion gives comparable ductility values for both. 

One important observation is that at a range of stresses above the Griffith criterion, 
where it is energetically favorable to fracture, the EP crack does not propagate, but the 
CLS crack does. This implies that lattice trapping from some source is important in the 
EP simulation. However, the stresses in the TB model (as computed by rigidly pulling 
apart two slabs) are comparable to those in the EP models. Therefore, something must be 
reducing the stresses needed to break a bond at the crack tip in the TB calculation. 

There are several possible sources for such weakening, including non-local effects 
and coupling between different strain components at the crack tip. 


WHAT IS CAUSING THE DIFFERENCE? 

Energetics; MO 

Rice criterion comparable 

Just, above Griffith criterion; 

EP crack, won't propagate 
TB crack will 

Stress induced lattice trapping; MAYBE 
Stresses in TB and EP comparable 

Even short-range TB (with higher decohesion stress) is brittle 
Hypothesis: 

decohesion stresses at crack tip are weakened in TB calculations 
Several possibilities; non-local effects, coupling 


Figure 18 
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WEAKENING: NON-LOCAL EFFECTS 


The Rice criterion (and its stress analog) are local. They assume that the relevant 
energies can be measured by creating the entire surface at once. However, in the real 
crack, a half surface is already present. It is possible that the broken bonds to the right of 
the crack tip weaken the bonds ahead of the crack tip. 

To measure this effect we created a penny crack in an otherwise ideal Si crystal, and 
rigidly pulled apart the two slabs. We monitored the forces on the atoms adjacent to the 
crack tip, and compared them with the forces on atoms in the bulk, far from the crack tip. 
In the EP calculations the forces on these two atoms were identical. Using TB the crack 
tip bond was indeed weakened. The difference in the peak forces between the crack -tip 
bond and the bulk bond was only about 2%, too small to explain the qualitative 
differences between the EP and CLS fracture simulation results. 

Other possible hypotheses for the reduced lattice trapping in the CLS simulations are 
currently being investigated. 


WEAKENING: NON-LOCAL EFFECTS 

Rice criterion (and stress analog): ‘ local'’ =► 

Assume entire surface is formed at the same time, but . . . 
In crack, half surface is already present 

Do adjacent, broken bonds weaken bonds ahead of crack tip? 



Create penny crack, pull apart, slabs, monitor forces on surface atoms 
EP: peak forces are identical 

TB: Force on crack-tip-bond atom weaker by < 2% 


Figure 19 
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CLS METHOD - FUTURE DIRECTIONS 


There are a number of limitations in the current implementation of the CLS approach 
coupling an empirical potential to a tight-binding description of a system. The first is that 
the code is limited to single component systems. The second is that Coulomb interactions 
due to the movement of electrons in the system are not included. 

We are planning a number of improvements to the code to address these issues. These 
include the extension to multicomponent systems, making it possible to study more 
complex, technologically relevant materials. Another is charge self-consistency, which 
will be even more important for heteroatomic systems where bonds are partially ionic. 
The last improvement is for full scalability of both the EP and TB region calculations on 
parallel computers, making it possible to study larger systems. 

The code is being developed and will be released through the DOD HPC CHSSI 
program. 


CLS METHOD - FUTURE DIRECTIONS 

Limitations In current implementation: 

• Single component systems 

• No charge self-consistency 

P I a n ned i m provernen ts: 

• Multicomponent systems 
more complex materials 

• Charge self-consistency 
more complex materials 

• Full scalability (EP. QM regions) 
larger systems 

CHSSI - DOD High Performance Computing software devel- 
opment project 


Figure 20 
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FUTURE APPUICATIONS 


Using some of these improvements, we are planning to use the CLS approach to 
study structural and mechanical properties. Planned applications include high strain-rate 
plastic deformation, as well as additional study of fracture. We hope to be able to 
simulate more complex, technologically relevant systems including other brittle 
insulators, e.g. semiconductors used for micro-electromechanical systems (MEMS) 
construction such as SiC and diamond, as well as high-toughness ceramic coatings. 
Although there are some technical challenges, we will also explore the application of this 
method to fracture in metals. Another field of application is friction and stiction, which 
are often the limiting factor for MEMS devices. Central to these phenomena is the 
interaction between surface chemistry and the mechanical loading. 


FUTURE APPLICATIONS 

Structural, mechanical properties; 

High strain rate plastic deformation 

Fracture in more complex, technologically relevant systems 

* Other brittle Insulators 

- MEMS semiconductors (Si, SIC., diamond) 

- high-toughness coatings: ceramics 

* Metals? 

Friction and stiction: (NSF NIRT - JHU, Naval Academy) 

♦ Big limitation for MEMS 

• Interaction between surface chemistry and mechanical loading 


Figure 21 
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CONCLUSIONS 


The concurrent coupling of length scales approach simulates a large system with a 
fast method with limited applicability, and where needed using a slower but more 
accurate method. The essential requirement for this approach is that the region where the 
slower method is needed be localized. 

The approach is most useful where there is no separation of time scales between the 
processes in the two regions. This leads to an interplay in dynamics between the two 
regions, making it impractical to compute the changing boundary conditions on the two 
regions a priori. 

The benefits of the approach are a combination of accuracy and speed. 


CONCLUSIONS 

Concurrent coupling of length scales: 

large system with a fast method, limited applicability 
where needed, use a slower, more accurate method 

Requirement: localized region for slower method 

Most useful when: 

• Not much separation of time scales 

♦ interplay In dynamics between regions 
boundary conditions on inner region changing, 
too complex to be computed ahead of time 

Benefits: combination of accuracy and speed 


Figure 22 
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OUTLINE 


This talk will start with a tutorial on computational methods and techniques. Two 
examples of atomistic calculations will be presented: mechanical properties and phase 
transformations. The author’s view of the future of computational materials science will 
be presented. 


THE FUTURE OF COMPUTATION 
MATERIALS SCIENCE: 

THE ROLE OF ATOMISTICS 

JAI 

'8 S' % Ssnt 

M. 1. BASKES 

LOS ALAMOS NATIONAL LABORATORY 

Methods 

Techniques 

Examples 

mechanical properties 
phase transformations 
Future 

/V 

N2SA 

• Los Alamos 

Figure 1 
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METHODS USED DEPEND ON THE ENTITY BEING MODELED 


The concept of entity is introduced. Each entity is connected to a specific method that 
may be used for calculations involving that entity. This talk with concentrate on semi- 
empirical methods. These methods have a basis in physics or chemistry, but use 
experimental data. 


METHODS USED DEPEND ON THE 
ENTITY BEING MODELED 

Method 

Entity 

Examples 

first principles 
semi-empirical 

electron 

density functional theory (DFT) 
quantum chemistry (QC) 

embedded atom method (EAM) 
modified EAM (MEAM) 

N-body, glue 

empirical 

, gram, 

dislocation 

Lennard -Jones (LJ), Potts 
dislocation dynamics (DD) 

phenomenological 


Ficks Law 
elasticity/plasticity 

jTt 

• Los Alamos 


Figure 2 
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ATOMISTIC MODELS HAVE TWO PURPOSES (I) 


There are two clear purposes for performing calculations at the atomic level. The 
first purpose is to obtain quantitative properties of specific materials for use in models 
using higher lever entities. Examples are given for first principles and semi-empirical 
methods. 


ATOMISTIC MODELS HAVE TWO PURPOSES ( 

• Obtain quantitative properties of specific materials 

- First principles 

- lattice constant 

- structural stability 

- elastic modulus 

- Semi-empirical potentials (fit to experimental data) 

- thermal expansion 

- melting point 

- yield strength 



• Los Alamos 


Figure 3 
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ATOMISTIC MODELS HAVE TWO PURPOSES (II) 


The second and equally important purpose of atomistic modeling is to obtain 
understanding of physical processes. Examples are given for empirical, semi-empirical, 
and first-principles methods. 


ATOMISTIC MODELS HAVE TWO PURPOSES ( 


Obtain understanding of physical processes 

- Model system (empirical) potentials 

- how specific properties affect collective behavior 

- dependence of yield strength on stacking fault 
energy 

- Semi-empirical potentials (fit to experimental data) 

- plasticity 

- phase transformations 

- First principles 

- diffusion 


Los Alamos 
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WE USED THREE WORKHORSE TECHNIQUES 


Techniques for atomistic calculations are presented. The requirements and output for 
each technique are summarized. 


WE USE THREE WORKHORSE TECHNIQUES 


«* Molecular statics 

- Requires total energy as a function of all atom positions 

- Convenient to also have forces 

- Yields minimum energy configuration 

• Molecular dynamics 

- Requires force on each atom as a function of all atom positions 

- Follows the atomic motion (vibrations) using Newton's second 
law 

• Monte Carlo 

- Requires total energy as a function of all atom 
positions 

- Yields a catalog of the most probable atomic configurations 


• Los Alamos 


Figure 5 
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WHAT MOLECULAR STATICS GIVES YOU 


The detailed information about what can be obtained from a molecular statics or 
energy minimization calculation is presented. 


WHAT MOLECULAR STATICS GIVES YOU 


• Structural Energies 

• Defect Formation Energies 

• Saddle Point Energies 

• Defect Geometries 

• Elastic Constants (0 K) 



• Los Alamos 


Figure 6 
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WHAT MOLECULAR DYNAMICS GIVES YOU 


The detailed information about what can be obtained from a molecular dynamics 
calculation is presented. Note that mechanisms and temperature dependent properties are 
obtainable. 


WHAT MOLECULAR DYNAMICS GIVES 
YOU 


* Mechanisms 

* Reaction Rates 

* Diffusion Rates 

* Interface Mobility 

* Elastic Constants (T) 

* Thermal Expansion 

* Thermodynamics 
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WHAT MONTE CARLO GIVES YOU 


The detailed information about what can be obtained from a Monte Carlo calculation 
is presented. Note that equilibrium properties are obtainable. 


WHAT IVIONTE CARLO GIVES YOU 


* Equilibrium at Temperature 

- positions 

- composition 

* Thermal Expansion 

* Thermodynamics 
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LIMITATIONS OF COMPUTATIONAL MATERIALS SCIENCE 


There are severe limitations in what computational materials science can do. Most 
important is our lack of information about underlying physical mechanisms. 


LIMITATIONS OF COMPUTATIONAL 
MATERIALS SCIENCE 

* Lack of understanding of underlying mechanisms 

- mechanical properties 

- phase transformations 

* Inability to query long enough times 

- processing - development of microstructure 

- properties - realistic strain rates 

* Inability to process large enough samples 

- boundary conditions 

- rnicrostructure 

* Inadequate potential interactions 

- multi-components 

• Los Alamos 

Figure 9 
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WE CAN USE ATOMISTIC CALCULATIONS TO HELP 
UNDERSTAND MECHANICAL BEHAVIOR 


We choose mechanical behavior as an example of atomistic calculations. The 
calculations will be used to investigate the strain gradient mechanism in plasticity. 


WE CAN USE ATOMISTIC CALCULATIONS TO 
HELP UNDERSTAND MECHANICAL 
BEHAVIOR 


• Classical materials science, physics, and mechanical 
engineering postulate that mechanical properties do 
not depend upon sample size 

• Recent experiments (torsion, indentation, bending) 
have shown that this postulate breaks down for 
small size samples 

• This behavior has been explained by postulating the 
importance of strain gradients 

• Atomistic calculations can investigate the “real” size 
mechanism 
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RECENT EXPERIMENTS SUGGEST A LENGTH SCALE IS 
INHERENT IN PLASTICITY 


Conventional theory = elasticity or plasticity theory. Experiments are torsion of small 
wires by Fleck and Hutchinson, indentation by, McElhaney, interfacial force microscopy 
(very fine scale indentation) by Houston et al. Large scale calculations were performed 
at Sandia. 


RECENT EXPERIMENTS SUGGEST A 
LENGTH SCALE IS INHERENT 
IN PLASTICITY 


• Conventional theory predicts that yield stress is 
independent of sample size (no inherent length 
scale) 

• Recent experiments show that yield stress 
depends on sample size 

• Molecular dynamics (MD) using the Embedded 
Atom Method will be used to investigate shear 
deformation in nickei 

• Calculations are performed using 100 to 100 
million atoms 
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GEOMETRY FOR ATOMISTIC CALCULATIONS 


The geometry we consider represents simple shear on a wire with rectangular cross 
section. 


GEOMETRY FOR ATOMISTIC 
CALCULATIONS 


* atoms follow atoms moved at constant velocity 


Newton’s 
second law 
temperature 
maintained at 
3O0K 



/ 


fixed atoms 


periodic, 
free surface, 
or GB 
in x 
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YIELD STRESSES BECOMES INDEPENDENT OF STRAIN 
RATE AT SMALL STRAIN RATES 


We show the yield stress normalized by the shear modulus vs. the strain rate. Note 
that the strain rates accessible in 111 MD are very large. The points represent the results of 
the MD calculations at various wire sizes indicated on the figure. The aspect ratio of the 
wire was held fixed at 2:1. The data shows that the yield stress becomes independent of 
strain rate for “low” strain rates. A simple inertial dislocation model is able to 
qualitatively describe this data. From the atom positions, we determine that the yield 
process is controlled by dislocation nucleation from the free (x-direction) surfaces. 


YIELD STRESS BECOMES INDEPENDENT OF 
STRAIN RATE AT SMALL STRAIN RATES 
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YIELD STRESS IS PRIMARILY DETERMINED BY SAMPLE SIZE 


Using the results of the MD simulations at the lowest strain rates, the shear stress, 
normalized by the shear modulus as resolved on the slip plane of maximum stress is 
plotted vs. a characteristic size, the volume to surface ratio. Also shown are the results of 
numerous experiments: interfacial force microscopy (IFM), indentation, and torsion. All 
of the data falls on an approximate straight line showing power law dependence of the 
yield stress on size. 


YIELD STRESS IS PRIMARILY 
DETERMINED BY SAMPLE SIZE 
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FACTORS THAT INFLUENCE THE YIELD STRENGTH IN FCC METALS 


The calculations presented above plus numerous other calculations lead us the 
following summary of factors that influence the yield strength in fee materials. Geometry 
effects, specifically size, are seen to be most important. Test conditions and material 
properties have a smaller effect. 


FACTORS THAT INFLUENCE THE YIELD 
STRENGTH IN FCC METALS 



Factor 

Extent 

Dependence 

Geometry 

Sample size 

1000 

-0.4 

X 


Aspect ratio 

10 

independent f > 4) 


Orientation 

4 

resolved stress 

Test 

Strain rate 

10 

independent f < s c (x) ) 

conditions 

Temperature 

1.2 

- 

Material 

Shear modulus 

10 

proportional to p 


Initial Dislocation 
density 

10 

“ 
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WHAT WE HAVE LEARNED ABOUT PLASTICITY MECHANISMS 


We may summarize what we have learned about plasticity from the MD calculations. 
In these small samples yield is determined by dislocation emission, not dislocation pileup 
as is commonly assumed. Deformation is found to be quite inhomogeneous even when a 
uniform strain is imposed. Finally we see a strong sample size dependence on the yield 
stress even in the absence of an externally applied strain gradient. These results force us 
to seriously question the strain gradient models in the literature. 


WHAT WE HAVE LEARNED ABOUT 
PLASTICITY fVIECHANISMS 


• Yield depends on sample size even in the absence of an 
applied strain gradient 

• Inhomogeneous deformation is the rule, not the exception 

• Strain rates In the inertial regime are important 

• The yield mechanism is dislocation emission fat least in 
small samples) - not dislocation pileup 
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WHAT WE HAVEN’T LEARNED ABOUT PLASTICITY 


Unfortunately MD atomistic calculations are not able to answer many questions of 
plasticity. These questions point out the clear need for modeling with different entities at 
large length scales. 


WHAT WE HAVEN’T LEARNED 
ABOUT PLASTICITY 

• What happens at low (realistic) strain rates 

- Other mechanisms, e.g. climb 

• Hardening 

- Restricted geometry (sample size) 

* Cell/microstructure formation 

- Limited sample size 

* Alloying additions 
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WE CAN USE ATOMISTIC CALCULATIONS TO 
HELP DESIGN NEW MATERIALS 


Our second example of atomistic calculations is in the area of material design. 


WE CAN USE ATOMISTIC CALCULATIONS 
TO HELP DESIGN NEW MATERIALS 


Pb-based solder has long been the workhorse of the 
electronics industry 

We now recognize the deleterious effects of Pb on 
the environment 

A major thrust has been initiated to replace Pb-based 
solder with a Sn alloy 

Atomistic calculations can predict properties of Sn 
alloys that are difficult to measure 
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GEOMETRY FOR REACTIVE WETTING CALCULATIONS 


The left figure shows a side view of a Cu substrate with a periodic array of small 
liquid Sn drops impinging on the surface. The right figure shows a top view of the 
sample. Modified Embedded Atom Method (MEAM) potentials are used to describe the 
interactions at the atomic level. These potentials describe the properties of Cu, Sn, and 
Cu/Sn alloys quite well. 



Figure 19 
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TIN REACTS WITH THE Cu SURFACE AT HIGH TEMPERATURE 


Using MD calculations, we see that at 700 K the Sn forms an ordered overlayer on 
the Cu (100) surface. At 800 K the Sn reacts with the Cu leading to a mixed layer. 
Experiment shows reaction at much lower temperatures. 


TIN REACTS WITH THE Cu SURFACE AT 
HIGH TEMPERATURE 
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A STEPPED SURFACE FACILITATES REACTION 


We investigate the effect of surface defects on the reaction temperature. The figure at 
the top left shows the initial defective surface. We now see that at temperatures as low at 
room temperature, the Sn reacts with the Cu substrate. 


A STEPPED SURFACE FACILITATES REACTION 
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WHAT WE LEARNED ABOUT WETTING MECHANISMS 


It is clear that surface defects are critical in the examination of the reactive wetting of 
Sn on Cu... It is satisfying to note that the MEAM potentials are able to represent the 
complicated processes occurring at the Cu/Sn interface. 


WHAT WE LEARNED ABOUT WETTING 
MECHANISMS 


* Liquid tin reacts with a smooth Cu (100) at 
temperatures greater than 800 K 

• Below 800 K an ordered layer of tin forms 


• A roughened surface reacts with liquid tin down to 
temperatures as low as room temperature 



Figure 22 
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WHAT WE HAVEN’T LEARNED ABOUT WETTING 


As in the previous case of mechanical behavior, it is seen that additional calculations 
are needed if we are to seriously approach the alloy design question. 


WHAT WE HAVEN’T LEARNED 
ABOUT WETTING 


• Solder alloy and surface impurity effects 

- Need potentials for Cu/Sn/X 

• Diffusion mediated mechanisms 

- Time limitations 

• Bond strength 



• Los Alamos 


Figure 23 
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DEVELOP MODELS BASED ON ENTITIES 
WITH INFORMATION PASSED FROM LOWER LEVEL 


We suggest that the most likely way for success in multi-scale modeling is the 
hierarchical scheme shown at the right. Calculations of the various entities are performed 
and information is extracted to be passed to the next level of computation. It is important 
to closely compare results to experiment at all levels of calculation. 


DEVELOP MODELS BASED ON ENTITIES 
WITH INFORMATION PASSED 
FROM LOWER LEVEL 


, dislocations 

electrons 

atoms 

Field 

% J| 

grains 

no 


add experimental 
information at 
every level 

retain only the 
minimal 
amount of 
information 

» Los Alamos 

Figure 24 
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IMPACT OF COMPUTATIONAL MATERIALS 
SCIENCE ON ALLOY DEVELOPMENT 


We feel that within the next two decades it is quite likely that computational materials 
science will be strongly influencing multi-component alloy development in the area of 
phase stability. In contrast the probability for strongly influencing alloy design through 
the prediction of properties is small. Somewhat intermediate is the ability of predict 
phase transformations and microstructure. 


IMPACT OF COMPUTATIONAL 
MATERIALS SCIENCE 
ON ALLOY DEVELOPMENT 


phase transformations 
microstructure evolution 

deformation 
fracture 

In order of increasing difficulty and decreasing 
probability of success 


Thermo-mechanical 

Processing 



» Los Alamos 
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KEY ISSUES FOR THE FUTURE 


We summarize the key issues for the future. 


KEY ISSUES FOR THE FUTURE 

• Reliable potentials - tradeoff 

- Reliability 

* Material classes 

• Multi-component systems 

- Computational speed 

• Information extraction 

- Connection to larger size (longer time) 
models 

- Visualization of mechanisms 


$ 


Changing entities - don’t do it ail with atoms 



» Los Alamos 


Figure 26 
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Materials Modeling and Simulation as a Technology Platform 
for Advanced Materials Development 


Gregg Caldwell 

3M Materials Modeling, Advanced Materials Tech Center 
St. Paul, MN 
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OUTLINE 


Outline 

• Materials Modeling Technology at 3M 

• Bottom Up 

» Atomistic Simulations of Organic Fluoromaterials 

* Predicting Solubility Parameters of Molecules. 

• Top Down 

* Hie Need for Property Predictions Via Molecular Simulations: 
Pavement Markings Performance. 

• The Middle Out 

* Self Assembly of Materials and the Mesoscale. 

• Summary. 

Materials Modeling 
Advanced Materials Tech Center 


Figure 1 
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MATERIALS MODELING TECHNOLOGY AT 3M 


The AMTC has significant expertise in modeling material behavior. The group is 
focused on the prediction of performance by using an integrated approach for various 
length and times scales - atomistic, molecular, polymer, mesoscale, fluid transport 
(including microfluidics), thru applied mechanics. The goal is to accelerate materials 
development with modeling and simulation and help develop new products and processes 
through computational methods. 


Materials Modeling Technology 

at 3M 

• Advanced Materials Technology Center 

* drives strong 3M global growth short-term and long-term 
through building and applying expertise in chosen materials- 
related technology areas, anticipating new technology platforms 
aligned with major 3M growth opportunities. 

• Materials Modeling Technology Group 

* accelerate 3M innovation by using mathematical and 
computational models of materials. 

* provide an integrated approach among the various length scales 
and time scales of materials behavior. 


Materials Modeling 
Advanced Materials Tech Center 


Figure 2 
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MATERIALS MODELING 


The modeling group within the AMTC use a hierarchal approach to materials 
modeling. The challenge of bridging the gaps in size and time scales is recognized with a 
“bottom up” and “top down” approach used to link well understood design and 
constitutive modeling approaches with first principles approaches. 



Figure 3 
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ATOMISTIC SIMULATION OF ORGANIC FLUOROMATERIALS 


Polytetrafluoroethylene (PTFE) is an important industrial polymer due to some of its 
remarkable properties. Other fluorinated materials also have significant industrial 
interest and the ability to predict structure and possibly performance of these materials 
would be of great value. 

Unlike the zig-zag shape of the polyethylene backbone, X-ray crystallography studies 
show 1 that polytetrafluoroethylene (PTFE) has a helical conformation due, in part, to 
larger van der Waals interactions among fluorine atoms. The first order phase transitions 
(20°C) below the melting point are believed to be due to a chirality transformation from a 
well ordered helical form to disordered ones and disordering among polymer chains. 

In order to look at this transformation a molecular dynamics (MD) approach is 
necessary. The accuracy of the MD calculations is only as good as the underlying force 
field. 

1. Bunn C.W. and Howells E. R. Nature, 174, 549 (1954). 

2. Mayo S.L., Olafson B.D., Goddard W.A., J. Phys. Chem., 94, 8897 (1990) 


Atomistic Simulations of 
Organic Fluoromaterials 

* Selected, a general force field for organic molecules: 

• Dreiding 2 . 

* The basic force field had limitations: 

• Did not deal with fluorinated systems correctly. 

* A collaboration between SWTs Materials Modeling Group 
and the Materials Process and Simulation Center led by Bill 
Goddard at the California Institute of Technology. 

Materials Modeling 
Advanced Materials Tech Center 


Figure 4 
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DREIDING EXP 6 FORCE FIELD 


Early experiments 1 demonstrated that the backbone dihedral angles in PTFE are 
displaced by 17° from a truly trans state (that is 0°) This can be attributed to the large 
non-bond interactions among fluorine atoms attached to the carbon backbone. In order to 
accommodate these unfavorable steric interactions, the carbon backbone must deform 
slightly. As a result, PTFE forms two stereo helical conformations. 

Full quantum mechanical calculations for PTFE, however, still remains a challenge 
due to the large number of atoms involved. However, molecular mechanics/dynamics 
methodologies have been widely used to model macromolecules and large systems. In 
order to carry out correct analysis from classical molecular methods, the critical issue is 
to obtain a highly accurate classical force field (FF) that models the atomic interactions 
among different atoms. 

As a first step, the torsional parameter of the Dreiding 2 force field was modified 
based on ab initio studies. 



Figure 5 
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QUANTUM MECHANICS 


To parameterize the torsional potential, we carried out a series of cib initio 
calculations for fluorinated alkanes from C 2 F 6 up to C 10 F 22 . (B3LYP with 6-3 1G* basis 
set). 

To capture the conformational characteristics classically, we used the standard 
Dreiding Force Field to describe covalent bonds, angles, and non-bond interactions. The 
van der Waals parameters for fluorine were taken from a previous study' 1 . The torsional 
potential was tuned to best represent helical displacements for fluorinated systems. The 
torsional potential used was a 6-term dihedral. 

3. KarasawaN, Goddard W.A , Macromolecules, 25, 7268 (1992). 


Quantum Mechanics* 


:w optimized force field. 
Based on. Dreiding 
Timed the torsional potential 
• 6 term dihedral in the form: 


^=ZZ^l(l + ^cos^,) 

■=1 i z 


Energy vs Dihedral Angle 



Dihedral angle 


Blue molecule is optimized 
using DFT. 



and Grey Molecule is 
optimized using new FF. 



Materials Modeling 
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VALIDATION - C20F42 


Perfluoro-n-eicosane (C 20 F 42 ) was used to validate the modified force field. Like 
other chain molecules, the crystal structure is composed of stacking layers. Within each 
layer, the molecules form a two-dimensional lattice. The phase transitions are first order 
and mainly due to dis-order in rotational motions within each molecule or translational 
motion among molecules. Perfluoro-n-eicosane has three solid phases 4 , namely, M 
(monoclinic) for T < 146 K; I (intermediate) for 146 < T < 200 K, and R (rhombohedral) 
for T > 200 K. 

4. Schwickert, H. J., Strobl, G., Kimmig, M., Chem. Phys. 1991 , 95, 2800. 


Validation - C 20 F 42 


• Crystal Structure: 

- At least 2 stable crystal phases 4 . 


H. Schwickert, G. Stroble, M. Kimmig, 
J. Chem. Phys., 95 (4), 2800 (1991) 


® A phase stable at ambient temperature. 

- rhombohedral phase (R) 

» Each unit cell occupied by 3 molecules 
» a :::: b =* 5.7 A, c =:: 85.0 A 
» a = 8 = 90°; y - 120 0 
» p = 2.16g/cm 3 

* A phase stable at low temperature. 

- Monoclinic phase 

» Each unit cell occupied by 2 molecules. 
» A = 9.65 A, b - 5.7 A. c = 28.3 A 
» (x = y = 90°; p ~ 97 0 
» p = 2.16g/cm 3 



Materials Modeling 
Advanced Materials Tech Center 
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VALIDATION - C 2 oF 42 LOW TEMPERATURE PHASE 


The unit cell in the simulation of the M-phase C20F42 consists of two molecules. One 
is left-handed, and the other is right-handed. After constant pressure energy 
minimization from molecular mechanics, the computed cell parameters are listed together 
with experimental data. 


Validation - C 20 F 42 


Low Temperature Phase 

Struct ure at minimization: 

• lxlxl : 

a=88.9°, (3=100.3°, y=94.6° 
a=9.73, b=5.63, c=28.71 A 
E=164.053 kcal/mol 
p = 2.24g/cm 3 

» Experimental 

» A = 9.65 A, b = 5.7 A. e = 28.3 A 
» a = y = 90°; (3 ~ 97 0 
» p= 2.23 g/cm 3 



Materials Modeling 
Advanced Materials Tech Center 
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VALIDATION - C 20 F 42 HIGH TEMPERATURE PHASE 


The high temperature phase having three molecules in the unit cell is much more 
difficult. Three different cases were tries: 

Case 1 : All the molecules in the unit cell were left handed. 

Case 2: All the molecules in the unit cell were right handed. 

Case 3: The molecules were alternating left and right handed. 

The results of the third case were the closet match to the experimental crystal 
structure after minimization. 


Validation - C 20 F 42 

high temperature phase 


• This is more difficult! 

.... How do the 3 molecules arrange 
themselves in the unit cell? 

- Crystal structure implies the 
ordering of the helices has 
disappeared. 

---• What is happening? 

Experimental: 

;i b - 5.7 A, c - 85.0 A 
<x- p - 90 c '; y “ 120 ° 

P 2.16 g/om 3 


I ried 3 Cases: 

•case I - all left-handed helices: 

ot=90.7°, p=89.3°, y=l 18.3° 
a=5.62A. b=5.60A, c=84.73A 

• case 2 - all right-handed helices: 

a=89.7°, p=90.3°, y=120.0° 

a=5.63A. b=5.72 A, c=84.8lA 

• ease 3 - alternative right-left helices: 
a=89.9°, p=89.9°, y=l 19.2° 
a=5.6lA, b=5.69A, 
c=169.8 A (equiv. C=84.9 A) 



(minimization) 
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VALIDATION - C20F42 


A 3X3 supercell of the M-phase was constructed (each layer is colored differently for 
clarity). The molecules are packed in a very ordered array. 




Validation — C20F42 

phase-3x3x3 


Looking down along Z-axis 


Looking down along Z’-axis 


zz 


Materials Modeling 
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VALIDATION - C20F42 


Molecular dynamics was then run at 120K (below the first transition) and at 240K 
(above the second transition). After 70ps of dynamics at the lower temperature, the well 
ordered crystal structure is maintained. When the temperature is increased, however, the 
crystal structure tends to disorder. The is due to a change in the helicity of each 
molecule. The barrier to change from a left handed helix to a right handed helix is ~ 
lkcal mol' 1 . During the higher temperature dynamics a bond changes form +15° to -15°. 
This in turn cases the other bonds in the molecule to switch thus changing its handedness. 
The adjacent molecular then changes handedness and the ordered of the entire crystal is 
changed. 


V ahdation - C 20 r 42 

M-phase-3x3x3 MD at 1.20K and 240K 


After 70ps of NPT 
dynamics at 120K 

(Looking down along Z-axis) 



X 


In the time average, the surface 
structure is maintained. 


After 70ps of NPT 
dynamics at 240K 

(Looking down along Z-axis) 



In the time average, the surface structure 
of the helices is smeared out => there is 
no ordering with regard to the helix 
handedness. 

^ ^ atef *ais Modeling 

Advanced Materials Tech Center 
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X-RAY DIFFRACTION PROFILE OF THE STRUCTURES FROM MD 


Simulated X-ray patterns averaged over the dynamics time period shows the 
difference in order in the molecules helicity at different temperatures. At low 
temperature (120K) the molecule has a lot of order in the higher 26 region. At higher 
temperature (240K) this order disappears as the molecules switch helicities. 



Figure 12 
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VALIDATION 


A recent paper' described the SAM formation of C 20 F 42 by vapor depositing the 
compound onto a glass slide. Meridial X-ray patterns of our simulated high temperature 
phase are in excellent agreement with the experimental X-ray patterns. 

5. Nishino T., Meuro M., Nakamae K., Matsushita M., and Ueda Y., Langmuir, 1999, 
15, 4321 


Validation. 



C 20 F 42 SAMS 

* C 20 F 42 vapor deposited onto glass 
(epitaxially grown ''single-like” 
crystallites) 


(0 0 3} 



3rd _ 24*. 

meridial 

reflections 
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PREDICTING SOLUBILITY PARAMETERS OF MOLECULES 


It is of great value to be able to estimate from first principles the Hildebrand 
solubility parameter for any given solvent or polymer. Predictions of this type may be of 
practical use for selection of polymers in blending, for understanding polymer kinetics 
and monomer distribution compositions in copolymers, or for the proper selection of 
solvents in new formulations and synthesis processes. 


Predicting Solubility Parameters 
of Molecules 

* in the literature, there sometimes as many values of 8 as 
there are researchers measuring them. 

* Distribution can be pretty large (especially for polymeric 
materials), 

* general rule of thumb - compare values from the same researchers only. 

* Can modeling give a quick, reproducible, way of calculating 
cohesive energy densities (CED) and solubility parameters? 

* A collaboration between the SWTs Materials .Modeling 
Group and the Materials Process and Simulation Center led 
by Bill Goddard at the California Institute of Technology 

Materials Modeling 

Advanced Materials Tech Center 
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SIMULATION DETAILS 


A simple hierarchal approach was again used in the simulation procedure. First, 
accurate charges were obtained for each molecule using quantum mechanics. Molecular 
dynamics were then carried out on an ensemble of the molecules and results averaged. 


Simulation Details 

® Atomic charges calculated via Hartree-Fock studies 
(6-3 1G** basis set) employing Jaguar. 

•••• Both Mulliken and Electrostatic Derived Atomic Charges Employed 
and Compared. 

• Molecular dynamics simulations carried out with the Cerius 2 
employing the Dreiding Force Field. 

• Molecular Dynamics Simulation protocol and 
CED/Solubility Parameter analysis carried out via the 
program CED developed Mario Blanco, Bill Goddard, and 
co-workers at the California Institute of Technology 

Materials Modeling 

Advanced Materials Tech Center 
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SOLUBILITY PARAMETERS 


A search of the literature found 39 compounds with reliable solubility parameters 
determined experimentally. The experimental points are from a collection of literature 
sources and there is quite a range of experimental values for each compound. 

Computational approaches to calculating the solubility parameter provide a consistent 
approach allowing much easier comparisons between compounds. 


Solubility Parameters 

9/ 


MuIIifeen 


Calculated and Experimental Solubility Parameters 
(39 Molecules) 



Average Experimental Value 


1,04 (cal/cm 3 ) 172 


| Materials Modeling 
* Advanced Materials Tech Center 
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PREDICTING SOLUBILITY PARAMETERS SUMMARY 


The method shows good agreement with experimental results and 3M is utilizing this 
as way of overcoming the need for time-consuming experiments. 


Predicting Solubility Parameters 
Summary 

* With Mulliken Charges, the average absolute agreement with 
experiment for 39 of the molecules is found to be 

1 ,04(cal/cm 3 ) 1/2 . 

• With Electrostatic Charges, the average absolute agreement 
with experiment for 16 of the molecules is found to be 0.88 
(cal/cm 3 ) 1/2 . 

* Simulations in progress to complete Cerius2/electrostatic 
charge simulations for comparison with Mulliken charge- 
based simulations. 

® Simulations in progress with Discover and the Compass force 
field to compare accuracy with the Dreiding Force Field. 

• Methodology is being incorporated into the research process 
at 3M to replace the need for time-consuming experiments. 

Materials Modeling 

Advanced Materials Tech Center 


Figure 17 
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PREDICTING SOLUBILITY PARAMETERS APPLICATIONS 


From an industrial stand point, predicting solubility parameters has many 
applications. 


Predicting Solubility Parameters 
Applications 

• Polymer/polymer miscibility for design of polymer blends 

• Polymer/polymer immiscibility for design of phase 
separating mat eri al s. 

• Polymer/ additive miscibility for optimization of additive 
selection. 

• Polymer/ solvent or moiecu le/ solvent miscibility for 
solvent selection. 

• Solvent/sol vent miscibility to design of solvent blends. 

• Molecular descriptor for development of Quantitative 
Structure Activity Relations. 


Materials Modeling 
Advanced Materials Tech Center 


Figure 18 


112 





NEED FOR PROPERTY PREDICTIONS VIA MOLECULAR SIMULATIONS 
PAVEMENT MARKINGS PERFORMANCE 


The Traffic Control Material Division of 3M Company is a recognized market leader 
in products for traffic management enhancing safety for all of us. This leadership comes 
from delivering innovative solutions to cities, counties, states and countries that need 
them. One of these solutions, pavement markings, utilized 3M retro reflective 
technology for lane markings. The customers want markings with road presence that 
means, durable markings with good adhesion and lasting retro reflectivity. 


The Need for Property Predictions Via 
Molec ular S imu I ation § 
Pavement Markings Performance 


* Road Presence 

- Adhesion 

---• Wear Resistance 
® Ease of Application 

- Wide Application Window 

- Minimum Road Surface 
Preparation 


* Visibility and Recognition 
- Color 

■ Retroreflective Performance 
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MATERIAL PROPERTIES FOR MARKING PERFORMANCE 


A major goal in the development of pavement marking tapes is conformance: the 
ability of the tape to maximize road contact and hence road adhesion. Intrinsic properties 
of the materials are needed as input ot continuum models based on finite element 
technology. The most important ones for our application are: Young's modulus, yield 
stress, glass transition, thermal and hygroscopic expansion coefficients, viscosity, contact 
angle with various pavement types and permeability both to water and oil. 

Our models have been used to predict the conformance of a particular tape in various 
roads. This slide has an example of the performance on two different roads with very 
different aggregates. Notice how performance depends on the type of road, conformance 
goes from 77% in one to only 33% in another one. 


Material Properties for Marking 
Performance 


Modulus Expansion Coefficient Tan Delta 

Yield Thermal Conductivity Contact angle 

Tg Viscosity Permeability 
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ACTUALS VS. PREDICTED ROAD CONTACT 


The chart summarizes results obtained with the finite element models. It compares 
predicted values with measured values in two different sections of the same road: the 
wheel track zone and the oil drip zone. The wheel track zone corresponds to the area of 
the lane where the wheels pass. The oil drip zone corresponds to the center of the lane, i. 
e, the area between the two wheel tracks. 

Our models capture correctly the ability of tapes to conform to various roads. 



Figure 21 
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MATERIALS MODELING 


New technology trends such as materials that interact, exchange or sense 
environmental changes, nanotechnology, microfluidics and micro machines are requiring 
information previously unavailable from traditional materials modeling methods, either 
fully atomistic or from a finite element side. 

To address this gap, mesoscale models with a characteristic size located between the 
atomistic scale and the continuum model scale are being investigated for numerous 
applications. 



Figure 22 
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SELF ASSEMBLY 


Of particular interest in soft-condensed materials is the self assembly of various 
components to form a structured fluid. The length and time scales at which these 
structures are formed are prohibitive for atomistic simulations, and mesoscale models are 
required to capture the rapid molecular kinetics as well as the slow thermodynamic 
relaxations of macro scale properties. 

By coarse graining the description of the molecules we gain orders of magnitude in 
both length and time-scales. 


Self Assembly 

V 






Seif-a.sse.mbly of fluids has length and time scales which are 


/e for fully atomistic studies. 


What is needed are tools that capture the the chemistry of the 
materials whilst coarse-graini eg the model. 


Atomistic 



Length: Nm 

Units: Atoms 

Time: Ns 

Dynamics: F-ma 


Mesoscopic 



100’s of nm (or more) 

Beads representing many atoms 
as much as milli-seconds 

Diffusion, hydrodynamics 
^ ^ ater * a ^ s Modeling 

Advanced Materials Tech Center 
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EXAMPLES OF MOLECULAR SELF ASSEMBLY 


Molecular self assembly is engineered for the development of a large variety of 
materials. Surfactants present themselves as unique building blocks. They are used to 
coat surfaces with monolayers. In solutions they aggregate into very complex structures 
that have been used as templates for other materials (such as porous systems) or as nano- 
sized packets (such as vesicles) or to catalyze specific reactions (such as carriers to 
enzymes). Molecular self assembly has also been used to carry out engineered nucleation 
for the formation of zeolites and clathrates. Further, ordering within zeolites pores or 
clathrate complexes can be used for separation or selective reactions. 


Examples of Molecular Self Assembly 




Self assembly on surfaces 


Nucleation & 
crystal growth 




VW “ 

n 
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1 
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formation of zeolites & 
clathrates 
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SELF ASSEMBLY IN MICELLES 


Specifically amphiles (or surfactants) provide a rich set of highly complex self 
assembled aggregates at various solution conditions. As one increase the surfactant 
concentration the monomers aggregate into spherical micelles. Increasing the surfactant 
concentration further, or adding electrolytes or various additives result in worm like 
micelles, bilayers, vesicles or bicontinuous phases. 


Self Assembly in Micelles 

qJ 



Concentration 



Figure 25 
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MESOSCOPIC SIMULATION - BLOCK COPOLYMER SYSTEMS 


Self assembling block copolymer materials have important commercial applications 
as thermoplastic elastomers, and compatibilizers in polymer blends. Recently more 
advanced applications for these materials are being developed. These include using block 
copolymers morphologies as templates for the production of nanostructures and in 
nanoscale lithography for example. 

While most blends of homopolymers are not compatible and show phase separation, 
block copolymers can only phase separate on a microscopic scale due to the connectivity 
of blocks. Depending on composition a large variety of micro-phase separated structures 
may be obtained. Mechanical properties of copolymers like the tensile strength may be 
significantly enhanced as compared to the corresponding homopolymers. 



Mesoscopic Simulations 

Block Copolymer Systems 

Block Copolymer are 

important industrial #< a dP# m. 

po,ymers - JfwS %% 

Used as compatibilizers. * 4# g| $ 8 j|i| fulpj 

- go lo the interface and | ^ \-'f 111 dt ' mr 

"blend” different ||f 1111® 

materials. ” 

Different blocks do not increasing A fraction 

“like” each other 


increasing S fraction 


separation. 
morphology drives 

• Final properties 

• Materials 
performance 

• Materials lifetime 
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THE PHASE DIAGRAM 


Like surfactants, changing concentrations of the various components of block 
copolymers, result in a number of possible morphologies. 

The polymers are represented by an ideal Gaussian chain model. 6 The non-ideality of 
the polymer is accounted for in an external potential. 

The free energy of this system can be obtained in terms of these external potentials, 
which in turn can be related to the density fields via a density functional relationship. 

In our group, we have been able to obtain computationally the phase diagram of a 
diblock copolymer. For this, we have used MesoDyn, a free energy method for blends 
and copolymer systems developed by Fraijee and implemented by Accelrys. Notice the 
right prediction of the various phases shown in the previous slide. 

6. www.accelrys.com 


The Phase Diagram. 
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TIME EVOLUTION FOR LAMELLAE FORMATION 


Typical processing times are orders of magnitude shorter than relaxation times. This 
results in non-equilibrium morphologies being trapped in the process. The systems are 
chemically the same, but have different material properties because various morphologies 
have been frozen in by processing. 

Being able to predict the morphology at various process times and the corresponding 
properties associated to that state would be a great benefit in industrial design. 

MesoDyn was also used to study the evolution of the A-B diblock copolymer toward 
its equilibrium morphology. 



Figure 28 


122 





SUMMARY 


Summary 

• Applied research studies are being earned out via atomistic 
simulations for: 

- polymer surfaces and polymer/polymer interfaces. 

- self-assembled molecular and polymeric materials. 

- organic/inorganic interfacial systems. 

- structure/property prediction. 



Hybrid atomistic/coarse grain methodologies are being 
incorporated into the 3M research program to expand the 
range of interfaces and materials that can be studied. 
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COMPUTATIONAL MODELING TO SUPPORT MATERIALS DESIGN 


The objectives of materials design may differ substantially from traditional objectives 
of life cycle engineering or assessment of durability. Rather than just selecting off-the- 
shelf materials, it is increasingly the case that we can consider alternatives among 
materials via computation. Since each new component and application has somewhat 
different functional requirements, it is desirable as well to develop the machinery that 
permits computational assessment of the relative suitability among a range of 
microstructures, or to be able to choose composition and process route that result in a 
microstructure that delivers the desired structure-property relations. This paper considers 
several projects at Georgia Tech that address multiscale modeling within the overall 
vision of materials design. 



Computational Modeling to Support 
Materials Design 


David L. McDowell 

Regents' Professor and Carter N. Paden, Jr. Distinguished 
Chair in Metals Processing 
Georgia Institute of Technology 
Atlanta, GA 30332-0405 


George W. Woodruff School of Mechanical Engineering 
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MULTISCALE MODELING- WHAT IS IT? 


In general, as one moves from increasingly fine length scales with high resolution, 
high fidelity mechanics models to more coarse scale representations at higher length 
scales, it is inevitable the (i) one loses description of dynamical state in favor of 
thermodynamical state and (ii) the number of model degrees of freedom decrease as one 
adopts lower order continuum models (compared, for example, to atomistic calculations). 
Note that although we normally speak of homogenization of properties over increasing 
length scales (cascade upwards), we must also understand the connection of macroscopic 
response to phenomena occurring at individual scales below this level, i.e., we must be 
able to deconstruct the multiscale description to extract salient statistics regarding 
evolving characteristics such as damage, inelastic deformation, and so on. 




Increase in DOF; move towards dynamical 
representation; computed or measured 


Georgia 

Tech 


George W. Woodruff School of Mechanical Engineering 


Figure 2 
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MULTISCALE MODELING- WHY DO IT? 


The most compelling argument to develop effective multiscale models is that relevant 
phenomena occur at disparate length and time scales, and are not easily connected within 
the umbrella of a single, unified theory. This is a real challenge in materials science and 
mechanics. Moreover we must account for microstructure variation if the ultimate 
objective is not just to derive first order average properties but rather properties that 
depend on extremal characteristics of microstructure. The modeling goes up and down 
through the length scales as we seek to explore effects of variation at different scales. If 
the goal is to design a fracture resistant microstructure, for example, we must be able to 
construct as well as deconstruct the hierarchical structure in order to isolate effects of a 
given scale. It is complicated by significant coupling of phenomena across length scales 
due to long range order in solids. We advocate a multiscale approach that respects 
conservation of mass, momentum and energy, and reflects similar characteristics of 
dissipation as we move back and forth between different types of models. The 
identification of schemes to capture evolutionary behavior over a representative volume 
in this regard is of immense practical importance and represents a fertile are for further 
development. The objective of a single, unified model for such a hierarchical, multiscale 
(in both length and time) is not terribly realistic. Rather, development of principles for 
linking models exercised at their most relevant and appropriate scale within some 
overarching set of heuristics or framework is the most desirable approach. 



!yf| j 1 I 


Phenomena occur at different length and temporal 
scales 

Fine scale behavior is often relevant to phenomena 
that depend on extremal characteristics - fracture, 
fatigue, turbulence; one-way homogenization 
unsuited for this task 

Practical limitations on computing time; DOF 
Limitations on confidence in local properties 


Georgia 

■ffech 


George W. Woodruff School of Mechanical Engineering 


Figure 3 
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MULTISCALE MODELING EXAMPLE - CAST ALLOY FATIGUE 


The first example illustrates the building of an hierarchical model for fatigue 
phenomena occurring over a wide range of length scales for cast A356-T6 alloys, ranging 
from several microns (individual Si particles) to 15-30 microns (fine scale gas pores) to 
30-90 microns (dendrite cells or secondary dendrite arm spacing) to shrinkage pores (up 
to 1 mm) and entrapped oxides tens of microns to mm range. 



Collaborators: 

M.F. Horstemeyer, Sandia-Livermore 

J. Fan, Alfred University 

K. Gall, Univ. Colorado 

Sponsors: 

Sub-contract with Sandia National Laboratories 
through CRADA with DOE/Big 3 
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Figure 4 
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REGIMES OF FATIGUE FOR CAST A356-T6 


It was found from elastic-plastic finite element analysis of individual cracked and 
debonded particles as well as ensembles of particles that the HCF-LCF transition 
corresponds to the percolation limit for cyclic microplasticity in the eutectic 
(interdendritic) channels. The transition is highly nonlinear. 



Regimes of Fatigue of 
Casting Alloys 


Si particles 
in A356-T6 




Limit 

Plasticity 


HCF-LCF demarcation is the 
percolation limit for cyclic 
microplasticity 
in the eutectic composition 
interdendritic channels 


j mktwlaaucity 


i mt 


Figure 5 
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FINITE ELEMENT SIMULATIONS 


Simulations were performed of debonded particles, including effects of particle- 
matrix contact during the compressive part of the cycle. This information was used to 
develop the structure of an R-ratio dependent correlation of cyclic microplastic shear 
strain over some characteristic length scale at the notch root (pore or particle). 



Particle 

Contact 


Debonded Inclusion;: 
3 Cycles, R = -1 ?: 


Local Maximum Plastic Shear Strain, y niax (%); 


Gall, K., Horstemeyer, M. F., McDowell, D. L., and 
Fan, J (2001), Int. J. Fracture, 108: 207-233. 
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Figure 6 
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PREDICTIONS OF THE OVERALL MULTISCALE FATIGUE MODEL 


Note two major aspects of the predictions for number of cycles to form a 1 mm 
fatigue crack. First, for various types of defects, the HCF life scatters over the observed 
experimental range for these types of defects. These calculations are deterministic for 
each type of defect - the scatter arises from variation of defect severity. Secondly, note 
the crossover in the LCF regime due to crack -particle and crack-pore coalescence. 
Higher levels of initial porosity are highly detrimental in the HCF regime, but are 
beneficial in the LCF regime. 





crack/particle coalescence 
j— crack/pore coalescence 


Hor&ertfatfy cad plate data 
1 2 evenly spacec fractured S> pastes 


McDowell, D.L., 
Gall, K., 

Horstemeyer, M.F 
and Fan, J., 
“Microstructure- 
Based Fatigue 
Modeling of Cast 
A356-T6 Alloy,” 
Engineering 
Fracture 
Mechanics, to 
appear, 2002. 


1 50 pore interacting wftfc 50 tor? pores 
and fractured $i parities 


Physics-Based Models Facilitate 

(a) consideration of competing effects 

(b) microstructure distribution effects 

(c) parametric studies - “what if’ games 

(d) optimization 


Figure 7 


133 





MULTISCALE MODELING EXAMPLE - CAST ALLOY FATIGUE 


With such models it is possible to identify clear physical mechanisms for various 
types of fatigue thresholds, shown below on the Kitagawa diagram. The three thresholds 
that arise in this model are (a) absence of cyclic microplasticity in the microstructure 
(yellow region), (b) propagation threshold for microstructurally small cracks based on 
CTD in A1 being less than one Burgers vector (below red line in the blue region), and (c) 
long crack propagation threshold (below blue line in the blue region. 



A356-T61 Al alloy 



Figure 8 
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MULTISCALE MODELING OF FRETTING FATIGUE OF TI-6AL-4V 


This work illustrates the use of microstructure-based models of a dual phase material 
in a multiresolution modeling scheme to better understand the phenomena of fretting 
fatigue. 



Chung-Hyun Goh, Georgia Tech 

J. Wallace & R.W. Neu, Georgia Tech (experiments) 

Sponsors: 

UDRI contract (admin: J. Burns and J.P. Gallagher) 
USAF Palace Knight program (R. Morrissey, co- 
advised by T. Nicholas) 

AFOSR (Monitor C. Hartley) 


or&f* = George W. Woodruff School of Mechanical Engineering 

Tech --- 

Figure 9 
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MICROSTRUCTURE FEATURES 


The specific Ti-6A1-4V microstructure to be modeled is a one with primary alpha 
(HCP) and lamellar grains with BCC beta phase and secondary alpha, shown below. 





Longitudinal 


A duplex microstructure: 

60% primary a grains and 40% 
grains with secondary a and p 
averaged in lamellar structure. 
a =930 MPa, a u =978 MPa 


Transverse 


[A duplex microstructure of Ti-6A1-4V] 
(Cortez eta/., 1999) 


George W. Woodruff School of Mechanical Engineering 


Figure 10 
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PLANAR TRIPLE SLIP APPROXIMATION FOR HCP ALPHA PHASE 


At room temperature, prismatic slip is the principal slip mode in the alpha Ti, which 
is modeled using a 2D planar triple slip approximation. Of course, the 2D approximation 
overconstrains specification of crystallographic texture. 




Figure 1 1 
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Specimen 


Specimen 


Applied stress, a 


COMPONENT FRETTING MODEL 


Using this crystal plasticity description, individual grains are modeled for the actual 
geometry of experimental specimen configurations, within the region of contact 
plasticity, while nonlinear kinematic hardening and elasticity laws are used in the 
intermediate and far fields from the contact, respectively. 


Figure 12 


Fixed 


Sponge ( 1 O' 5 /:) 


Mesh Size (f>) 1 

Radius of Cylinder (R) 5000 


Planes of Symmetry 
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COMPARISON TO FRETTING FATIGUE EXPERIMENTS 


The first example illustrates the building of an hierarchical model for fatigue 
phenomena occurring over a wide range of length scales for cast A356-T6 alloys, ranging 
from several microns (individual Si particles) to 15-30 microns (fine scale gas pores) to 
30-90 microns (dendrite cells or secondary dendrite arm spacing) to shrinkage pores (up 
to 1 mm) and entrapped oxides tens of microns to mm range. 




Georgia 

Tech 


George W. Woodruff School of Mechanical Engineering 


non-propagating 

cracks 


Similar... 


Figure 13 
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INTERFACE SEPARATION LAWS FROM MOLECULAR CALCULATIONS 


The objective of this work is to embed results of nonequilibrium processes (i.e., 
dislocation generation and motion) into internal state variable continuum models for 
cohesive interface separation potentials to treat interfaces with one or more ductile 
metallic constituent. Initial applications are to twist and tilt grain boundaries in Cu. 
Rather than using some kind of overlapping matching scheme, we use the molecular 
dynamics and statics results as “numerical experiments” from which we extract 
information regarding the form of the potential for combined normal and shear 
separation. To date, such potentials are purely phenomenological and typical have 
nonlinear elastic character. Some recent developments have included the notion of a 
degradation of unloading stiffness, motivated by secant stiffness behavior of brittle 
microcracked solids. Such approaches are likely not descriptive for interfaces between 
ductile metallic constituents. 



Interface Separation Laws from 
IB Simulations - Stress State and 
History Effects 


D. Spearot, K.l. Jacob and D. L. McDowell 
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Figure 14 
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COHESIVE FINITE ELEMENT METHOD - A MEANS FOR ANALYZING 
RANDOM MICROSTRUCTURES 


Distributed cohesive fracture interface elements have proven to be a useful tool for 
modeling the influence of microstructure morphology of multiphase brittle material 
systems. We cite below some aspects of the approach. Needleman (1987) was the first 
to use such approaches in computational schemes as an alternative to conventional 
fracture mechanics. The advantage of using such models is that distributions of cracks 
and crack growth paths are dictated naturally by microstructure without user intervention. 
Of course, there are issues related to convergence, mesh refinement, scaling from 
nanoscale separations to micron-scale separations, path history dependence, modeling 
cohesive versus adhesive fracture processes, modeling bulk (not just interface) fracture, 
and artificial compliance enhancement due to distributed separation elements through the 
microstructure. 



• Cohesive Finite Element Method (CFEM) 


• No assignment of fracture path 

• Consideration of significant, 
representative volumes 

• Facilitates the introduction of lower 
length scales into the analysis of fracture 


Individual Cohesive surfaces along 

constituent element boundaries 
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Figure 15 
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MOLECULAR DYNAMICS CALCULATIONS 


The embedded atom interatomic potential is used in dynamical calculations of 
separation of idealized Cu-Cu grain boundaries, as shown below. We have considered 
effects of both periodic and non-periodic boundary conditions. 





: 


Cu grain boundary interface 


45 degree tilt 
lattice 

misorientation 


Boundary conditions 

• Last plane in +/- [010] direction is 
frozen to remain planar 

• Tensile displacement rate of 1 A/ps 
applied ‘fixed’ planes 

• Temperature = 300 K 


George W. Woodruff School of Mechanical Engineering 


Figure 16 
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PATH HISTORY DEPENDENCE OF NORMAL AND SHEAR SEPARATION 


Significant differences are observed between normal and shear loading (deformation) 
paths for the interface. In particular, in tension with periodic unloading, the average 
stress returns to the origin. In shear, by contrast, there is a residual deformation of the 
representative volume of atoms, indicative of irreversible dislocation generation and 
motion. The degree of reversibility of the displacement in the normal separation depends 
on the scale of the characteristic volume of atoms considered - with relatively small 
volumes (thousands of atoms), the image forces on dislocations emitted from the 
interface are high enough to drive them back out of the surface as separation proceeds. 
This is not the case with shear. 



* Atomistically inspired damage- Atomistic calculations 
dependent separation show that normal and 

tangential responses to 

f f damage are different 

* * ' BA, 


C a. 4 S" 


Cu 4ti Latfc e« 



Figure 17 
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MD-INSPIRED CONTINUUM SEPARATION MODEL 


The elements of a continuum model with structure extracted from these kinds of 
molecular calculations are akin to the elements of a conventional plasticity model. A 
hyperelastic separation potential is still assumed with regard to a relaxed configuration of 
separation, but internal state variables are introduced that relate to the degree of interface 
damage within some boundary layer (based on changes of coordination number) and a 
measure of dislocation density on various slip systems (ascertained from centro- 
symmetry parameters - still in progress). The interface separation normalizing 
parameters are assumed to depend on these ISVs. The ISV evolution equations are based 
on the MD calculations, and are size scale and interface structure dependent. A challenge 
is the scaling of the continuum model to address modeling of fracture processes at 
various length scales. 



Embedding atomic characteristics into 
phenomenological model 
••• Effective nonlinear elastic separation potential 


■j is a set of interface attributes (ISVs) 


The interface topology and structure are embedded 


Associated evolution equations 




Figure 18 
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PRINCIPLES FOR RIGOROUS EQUIVALENCE OF CONTINUUM AND 
PARTICLE SYSTEMS FOR MULTISCALE MODELING 


One of the problems with existing propositions to link atomistic or other discrete 
particle system calculations with continuum fields is that many approaches seek to match 
only certain aspects of each solution. In this work, we propose a theoretical structure, 
based on the dynamic principle of virtual work, to provide transfer of all relevant 
information from the particle system to the continuum. In so doing, we effectively use 
the discrete particle system solution as the constitutive law, and construct current 
configuration continuum fields that are dynamically consistent. After doing this, one 
may then subject the resulting continuum fields to spatial and temporal averaging 
procedures, but these are relatively well established in the mechanics literature. It is a 
separate task to develop continuum constitutive equations. 



Continuum-Discrete Wolecular Dynamic 
System Equivalence 


M. Zhou and D, L. WcDowel 


M. Zhou and D. L. McDowell, Equivalent Continuum for Dynamically 
Deforming Atomistic Particle Systems, submitted to Phil. Mag. (2001), 
revised December 2001. 
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Figure 19 
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ELEMENTS OF DYNAMIC EQUIVALENCE OF CONTINUUM AND 

PARTICLE SYSTEMS 


In addition to many desirable attributes of the resulting continuum, work-conjugate 
stress and couple stress fields emerge (the latter only in the case of non-central 
interatomic moments). As usual, for example, the Cauchy stress is the work conjugate to 
the rate of deformation tensor. We also assert that if interatomic forces are purely of 
central type, then any couple stresses must arise from higher scale morphological 
arrangements of defects. The assessment of the equivalent continuum may be performed 
either in parallel with the atomistics or as a post-processing step, and is computationally 
rather intensive. None-the-less, the proposition provides a vehicle to speak of rigorously 
(rather than heuristically) equivalent fields. 



Elements of Fully Dynamic Equivalence 



* Equivalent continuum (EC) in fully dynamic sense: 

♦ Equivalence of internal, external and inertia work rates on 

arbitrary scale 

♦ Conservation of momenta, energy and mass 

♦ Definition of work-conjugate deformation field 

♦ Account of nonlocal effects 

♦ Explicit account of structure 

♦ Account of heterogeneity 

* Work-conjugate stress and couple stress fields: 

♦ Nonpolar materials: Couple stress unnecessary at nano-level; 

arises from scaling 

* Solution is computationally intensive 


Georgia 
© “Ifech 
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NONLOCAL EFFECTS OF ATOMIC INTERACTIONS 


Within some specified surface in the current configuration that passes through a series 
of atoms, interactions with atoms outside the surface are treated as body force and body 
moment fields. Hence, nonlocal action is dominant for small volumes of atoms, and is 
less significant as the ratio of the size of S to cut-off radius for lattice interactions 
increases. 
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“CONTINUUMIZATION OF FIELDS” 


The continuum is imbued with continuous mass density characteristics and 
continuously distributed stress and couple stress fields based on interpolation functions 
akin to the finite element method. 



Figure 22 
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WORK EQUIVALENCE OF CONTINUUM AND DISCRETE 
PARTICLE SYSTEMS 


Below are shown the dynamic work equivalence of internal forces (in black), external 
forces (in blue) and inertial forces (in red) for a current configuration volume element 
enclosing a finite set of particles/atoms and the equivalent continuum. It is noted that 
continuum stress, deformation and density fields are all expected to vary substantially 
over distances comparable to interatomic spacings, as they mimic the atomistic solutions. 
However, these continuum fields may then be subjected to classical volume and time 
averaging procedures, as done, for example, in turbulent flow. 



Particle System and Continuum Equivalence 
(On Arbitrary Size Scale) 


Continuum 


Discrete atomistic 


Gi wtfge W. Wiwidruff School <>f Kngtnemttg 


Tech 


Figure 23 
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Quantum Simulations of Matter: 
solids, liquids and nanostructures 


Ousntuas Simulations 


Stats of SSro Art Exfjerlments 




Figure 2 
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Computation as a tool for scientific 
investigation 


• Why ab-initio simulations? 

• Some interdiscipl in ary problems where ab-initio 
simulations can have an impact 

• How the methods developed in the solid state physics 
community have evolved and where they have ' spread' 

• How these methods ‘join’ quantum chemistry methods 
and how they are reaching out to the biology arena 

• Major limitations and challenges 


Figure 3 



Microscopic modeling based on classical molecular 
dynamics may have limited predictive power 


[ Ctessfca! 

• The materials properties are computed 

j Moteeutar Dynamics 

assuming a given interatomic potential 


fitted to experiment. 


• Fits are done using data taken over a limited 


range of thermodynamic conditions 


• No explicit description of bond formation and 

jp 

breaking -> no chemistry 

| taswton sqyjj&ohs w tom 

•No explicit description of electronic properties 


Complex states of matter involving forming 
and breaking of bonds and involving different 
phases cannot be described with classical 
molecular dynamics 



Figure 4 
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Ab-initio simulations require only 
fundamental physical constants as input 



Quantum simulations are a robust and sound tool for predicting 
properties of materials which cannot be measured directly 


Figure 5 



Ab-initio molecular dynamics is numerically 
intensive 


Classical MD 


Ab initio MD 


At each MD step: 

Compute forces from a 
given model energy 
function 

Update atomic positions 


At each MD step: 
Solve the Schrodinger 
equation for electrons 
Derive forces from 
electronic wavefunctions 
Update atomic positions 


Millions of atoms 


Hundreds of atoms 


Effective ab initio MD simulations require the use of supercomputers 



Figure 6 
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Ab-initio molecular dynamics simulations contain 
approximations 


Staoironic Structure 
fat gte-cttonc} 


fcSolecufcsr Oynamtes 
;pMa» Bqm tlen* fa* tern) 


•$• is the* mechanical 

giOMfHi stetfc emrgy of N 
tttlemcth'g oteehwrs in the 
Retd of tons 


Figure 7 



Quantum Monte Carlo methods provide a near exact \ 
solution of the Schroedinger equation 


Variational Monte Carlo 


Diffusion Monte Carlo 



Provides near exact solution to the many-electron problem 

Quantum Monte Carlo (QMC) methods are 
computationally very demanding 
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A new linear scaling QMC code has been 
developed based on Wannier functions 


Standard applications of QMC scale as N 3 where N is the 
number of particles 



0 200 400 600 800 1000 


Number of electrons 

The linear scaling of the new algorithm allows us to handle much 
larger systems (A.Williamson, R.Hood, J. Grossman, PRL 2001) 

Figure 9 

Size and time scales in ab-initio 
molecular dynamics 




•'Solving' size problem (N 3 scaling) is not enough for ab- 
initio MD : larger systems require in general longer 
simulation times. 



TB-O(N) with ~ 5000 atoms and ~ 100 pS Current major 
(A. Canning, J.Kim, G.G. PRL 97, JCP98). limitations of O(N) 

methods for MD: 

Limitations from slow motion of polymeric- 

, , -Convergence 

like fullerene structures 

-Same accuracy for 

T.Frauenheim, M. Sternberg, G.G.. Comp Phys. Com. 1999; different’ 

G.G. Phys.Stat.Solidi (2000); Comp. Mat. Ssci 1998. COnficjur3tions 


Figure 10 
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From solid state physics to chemistry and f(| 
biology' 


• 'Standard model of solid’ (Kohn-Sham, 1 965 -> M. Cohen et al. 1 979 ) 

• Combine statistical mechanical methods (Molecular Dynamics) and 
the 'Standard Model for Solids’ (Density Functional Theory): R.Car and 
M.Parrinello, 1985 -> 

• Improve Density Functionals (Gradient Corrections): ca. 1990 from 
condensed matter physics to chemistry 

• Quantum Monte Carlo calculations for real materials (D.Ceperley et 
al.) : ca. 1990 -> 

•Coupling classical/quantum methods towards biology. 


Figure 1 1 


Progress in condensed matter physics, materials 
science and chemistry using computations 



Progress in ‘traditional materials modeling’ at different 
length scales: predict new materials, predict new 
properties, “maneuver atom by atom” 


Progress at a “new” length 
scale, in the nanosdenes world 


microscopic materials 
modeling into the 
biology world 


Figure 12 
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Figure 14 
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A numerical fracture experiment 




Quantum mechanical simulations of 
micro-fracture in a complex material 


The computer generated 
a-SiC network is 
representative of a real 
material {stoichiometric CVD 
grown sample) 

We learned about: 

— Onset of microfracture 
and chemical properties 
— Elastic and hardness 
properties 

— ‘Definition’ of surfaces 


* <3. GstHI, F. Gygi> A. Catelteni. Phys. Rev, Lett, 1939. |, 


PC ) >4* v-6t -f 


Figure 16 
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Elastic and hardness properties of a-SiC 


Young modulus (fifa) 


First ab-ir»itio determination of hardness 
of a disordered semiconductor 


M. A. El Khakaot et at, J. Mat Res. 9, 96 (1894); Phys. Rev. B SI. 4903 (1985) 


Figure 18 
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We use a unique combination of quantum 
simulations tools to study nanostructures 



Quantum Monte Carlo : 
Optical Gaps 



Different 
accu racy 
for 

different 

properties 


1 # 

pip : 




:> •> 

, 

Many-Body Perturbation: 

Full Optical Spectra 



Ab-intio Molecular Dynamics: 
Structures, Dynamics, Phononsl 


Empirical Pseudopotentials: 
Quasi-Bulk Systems 



Figure 20 
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Ge nanoparticles : separate quantum 
confinement and surface reconstruction effects 



Dots with diamond and Tetragonal Structures 


We have studied the structure of the 
crystalline core and the surface 
reconstruction of Ge clusters up to 
2.1 nm (L.Pizzagalli, G.Galli, J.E.KIepeis 
abd F.Gygi. PRB2001) 


Figure 21 


Surface chemistry of Si nanoparticles 



We have studied the influence of 
different surface passivants (O, r 

S,F,CI,0-H and CH2) on the optical 
gap of hydrogenated Si 
nanoparticles, ranging in size from 
1 to about 2 nm (from Si35H36 to 
SM48H120). 



A. Puzder, A. Williamson, J. Grossman 
and G.Galli, PRL2002. 


„V 


'4 ■ x 




We have used a combination of Density 
Functional and Quantum Monte Carlo 
(QMC) calculations. 


QMC was used to verify LDA trends as a 
function of size, and to carry out a 
quantitative comparison with 
experiment. 


Figure 22 
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Figure 23 


Figure 24 


Surface chemistry of Si nanoparticles 


We found that a double bonded 
oxygen atom on the surface 
greatly reduces the single particle 
and optical gap for hydrogenated Si 
dots with diameters less than 2 nm. 

For sizes larger than 2 nm, oxygen 
states disappear from the gap and 
the HOMO-LUMO states are Si 
core states. 

These results provide a consistent 
interpretation of several 
experiments appeared in the recent 
literature (Van Buuren 1998, 
Schlupper 1994 and Wolkin 1999). 


Si8?H?s 


1.S nm 


A. Puzder, A. Williamson, J. Gross 
man and G.Galli 
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Complex systems can be modeled using ab- 
initio molecular dynamics 



Nanostructures 


Hydrogen under 
pressure 


Semiconductor 
turfaces (SiC) 
and quantum 
dots (Si, Ge, C, 


DNA model 
backbone in 
solution 


Water and 
hydrogen bonded 
water mixtures 


Structural, Electronic 


and Dynamical Properties 




Figure 26 
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How does DNA interact with its environment? 


J.White.E.Schwegler, G.Galli and F.Gygi J.Chem.Phys. (2000) 

E. Schwegler, G.Galli and F.Gygi Phys Rev. Lett. (2000). 

F. Lighstone, E.Schwegler, G.Galli and F.Gygi Chem.Phys.Lett (2001) 
E.Schwegler, G.Galli and F.Gygi Chem.Phys.Lett. (2001) 


Figure 28 
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Model of the DNA backbone 



Gauche-Trans 


Trans-Trans 


* Smallest realistic model of the 
phosphodiester linkage 


Contains hydrophobic, polar, and 
ionic solvation sites 




Figure 29 



The importance of counter ions on 
DNA backbone structure 


Contrary to long held assumptions 
based on gas phase data, Na* plays a 
crucial role In driving a conformation 
transition of dimethyl phosphate 


This transition can have important 
implications for DNA structural stability 


E.Schwegler, G.Galli and F.Gygi Chem.Phys.Lett. (2001) 


Figure 30 
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E.Schwegler et al 


Work in progress: understanding the 
reactivity of anti-cancer drugs 


Nitrogen mustard based DNA alkylating agents 


Activation steps involve intramolecular cyclization reactions 
that are poorly understood 


Figure 3 1 


Work in progress: understanding the reactivit| 
of anti-cancer drugs 


Figure 32 


Cyclization reaction in solution using constrained ab-initio 
molecular dynamic simulations 
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Grand challenges 


Time scales and size effects 


Identifying and 
modeling rare events 
{e. g< proton transfer, 
drug-DNA Interactions) 


Including electronic 
excitation properly 
(e. g. hydrogen) in a 
dynamical way 


Define proper analysis 
tools and take full 
advantage of electronic 
structure data 


Need to form scientists who know how to work 
at the interface of different disciplines 


Figure 33 



I would like to acknowledge my collaborators: 



Figure 34 
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OVERVIEW AND 

LIST OF COLLABORATORS AT VARIOUS STAGES OF THIS WORK 


Cleaning solutions, LCDs, and biophysical systems all rely on structured solutions for 
key aspects of their behavior. Simulations have a useful role to play in examining, 
learning about, making predictions concerning, and utilizing these systems. Atomistic 
simulations can be fairly routine but even if enormous resources are used in such studies 
many important phenomena he outside their range of applicability. As a result simplified 
models with dramatically higher computational efficiencies have been used for many 
such studies. Generally, such approaches have used generalized models that do not 
correspond well to specific systems. Here we describe largely successful attempts at 
creating models with mimic specific systems fairly well. 
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OUTLINE 
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THE FOCUS OF OUR STUDIES 


Our focus is currently on structured aqueous solutions such as those formed by 
surfactants or phospholipids. 

The types of phenomena of most interest to us generally involve collective behaviors 
that are beyond what can be studied using atomistic simulations. 


Class of Problems 

Structured solutions including 
solutions containing common 
surfactants 

Types of Phenomena 

- Phase Behavior 

- Solvation Behavior 

- Adsorption at Interfaces 


Figure 3 
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EXAMPLES OF STRUCTURED SOLUTIONS 


These are images from atomistic simulations of a commonly studied surfactant SDS 
(sodium dodecyl sulfate). The molecules in our studies consist of a hydrophilic group 
attached to one or two hydrophobic groups. The competing drives to maintain contact 
between the water and the hydrophilic groups while minimizing water-hydrophobic 
contacts leads to aggregation into structures with spacial extents limited by the size of he 
molecules involved. The top image is of a roughly spherical micelle while the bottom is 
from a lamellar liquid crystalline phase consisting of bilayers of surfactants. Colors: 
small yellow, sodium ions; large yellow, sulfur; red, oxygen; gray, hydrogens and 
carbons. Water is present by not shown in the micellar system. 



Figure 4 
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SIMULATION METHODOLOGIES 


Different levels of simulation methodology that may be applied to surfactant systems. 
Quantum mechanical studies, involve solving Schrodinger’s equation for the electrons in 
the system. Atomistic models generally treat atoms as spheres which may be linked 
together to form molecules. Electrostatic interactions are included by placing fractional 
point charges on the nuclei of the atoms. Coarse grain models (not shown) involve 
representing local groups of atoms that are fairly rigid by single sites. Such sites can also 
be linked together to represent molecules. Mesoscopic models are similar except that 
they generally represent many atoms with considerable internal flexibility. A number of 
continuum approaches exist, including computational fluid mechanics. The amount of 
computer resources needed to study a given volume of solution dramatically decreases as 
one goes from quantum, through atomistic, coarse grain, mesoscopic to continuum 
methods. 


Simulation Methodologies 



Figure 5 
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SIMULATION CAPABILITIES VS. LENGTH AND TIME-SCALES FOR SDS 


This is a rough plot of the regions of applicability in study duration and system size 
for various computational methods. As well, key time and length scales for SDS are 
indicated. Many important phenomena in these solutions are related to aggregate 
formation and evolution and currently require coarse grain or mesoscopic methods. 



Figure 6 
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ATOMISTIC SIMULATIONS 


Applications of atomistic simulations to structured solutions is illustrated in the 
earlier images for SDS and here for a commonly studied phospholipid, 
dimyristoylphosphatidylcholine (DMPC). These studies almost always require that the 
system be constructed with roughly the right structures prior to conducting the simulation 
itself The models employed for this type of system are fairly mature these days and 
many properties can be reliably estimated from such simulations. 



Figure 7 
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WHAT CAN BE LEARNED FROM ATOMISTIC STUDIES? 


While atomistic simulations can be used to study local structure and dynamics, length 
and time-scale limitations force us to consider other methods. 


What can be learned from 
atomistic studies? 

Local structure and dynamics can be well 
characterized: 

•chain conformation, 

•association of small species 
•diffusion of small species 
•surface tension/surface area per lipid 

Limitations 

•System size (<100 A) is smaller than 
many phenomena of interest. 

•Time simulated (<10 ns) is generally 
too short to study many collective 
behaviors. 

Cannot study large-scale motion and 
reorganization. 


Figure 8 
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COARSE GRAIN 


Coarse grain models are the next level of simplification above atomistic models. The 
use of coarse grain models in studying structured solutions is not new. While such 
studies generally involve simple models that do not mimic specific systems, they have 
shown that these models can capture many of the phenomena present in the real systems. 
Our work is an experiment to see if coarse grain models can be parameterized to 
accurately mimic specific systems (chosen in advance). To do this we force the models 
to reproduce key structural features from experiment and atomistic simulations. We will 
describe how such a model was constructed for DMPC. 


Coarse Grain 

To Get Beyond The Limits of Atomistic 
Models Simplify The Representation. 

Coarse grain models are the Next Level up in 
simplification. 

- General phenomena reproduced fairly well 
(Gary Larson) by simple models 

- No serious attempt to make such models 
accurate. 


- To see if simplified models can be made 
accurate. 


- Force model to reproduce key structural 
structural features from experiment and 
atomistic simulation (similar to A. 
Lyubartsev’s work). 


Figure 9 
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SOME ASPECTS OF THE MODEL ARE DERIVED FROM SEPARATE 

SIMULATIONS 


As much as possible a divide and conquer approach is used in parameterization of 
coarse grain models. Models for bulk water and hydrocarbons can be constructed 
separately. Each coarse grain water site represents 3 water molecules. Similarly 
hydrocarbon site represents 3 carbon atoms and their associated hydrogen atoms. 
Hydrocarbon sites can be linked together to represent hydrocarbon chains. The 
parameters used in the models were adjusted to reproduce the bulk densities and vapor 
pressures of pure water and hydrocarbons as determined by experiment. As well as the 
average length and stiffness of hydrocarbon chains as estimated by atomistic simulations 
is reproduced. 


Some aspects of the model are 
derived from separate simulations. 

The geometry of the hydrocarbons is 
determined by comparison with the results 
of atomistic simulations of bulk hydrocarbons. 

Hydrocarbon and water models are also 
based on the density and vapor pressure 
of the bulk liquids. 



Figure 10 
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WATER/HYDROCARBON INTERACTIONS 


The parameters governing the interaction between water and hydrocarbon sites were 
adjusted until the width of the water-hydrocarbon interfaces was about 5 A (roughly the 
value found in atomistic simulations). 


Water/hydrocarbon interactions adjusted to 
approximately match the width of the 
water/hydrocarbon interface. 



Figure 1 1 
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SITE MAPPING FOR THIS COARSE GRAIN MODEL 


The mapping atoms to sites for the coarse grain model is not unique. Generally, it 
involves identifying groups of atoms that are essentially rigid and representing them by 
single sites in the coarse grain model. 
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HOW TO MAKE THIS WORK? 


There are many parameters in our coarse grain model for DMPC, nearly all of which 
are determined by comparison with atomistic simulations. The procedure involves 
choosing an initial set of parameters, performing a simulation using the coarse grain 
model, comparing the results with a related atomistic simulation and adjusting the 
parameters in the coarse grain model to improve agreement. This procedure is repeated 
many times until acceptable agreement is obtained. The parameters controlling the 
internal geometry of these molecules converge quickly. Non-bonded interactions 
between sites are parameterized based upon the radial distributions functions, g(r) 
(essentially, the density of a type of site as a function of distance from the site of interest, 
normalized by dividing by the average density). 


How to make this work? 
Strategy: 

Try to match atomic level structure as 
closely as possible using the coarse grain 
model. 

Choose bond lengths (14) and angles (16) 
and their force constants appropriately. 

Non-bonded interactions: 

Electrostatics (2) and van der Waals 
potentials (28) chosen to match atomistic 
g(r)’s. 


Figure 13 
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INTERACTION POTENTIAL 


Most g(r)’s are relatively unstructured and so the corresponding potential can have a 
simple form (e.g. a Lennard- Jones 9-6). However, g(r)’s amongst the hydrophilic groups 
are highly structured in a complex manner. For these interactions, we use tabulated 
potential (following Lyubartsev’s approach). These potentials are adjusted until the 
corresponding g(r)’s match closely. In the top figure the atomistic (solid lines) and 
coarse grain (dashed) g(r)’s for the choline group are given. The bottom figure gives the 
corresponding tabulated potentials obtained from the iterative fitting procedure. 




Distance in Angstroms 


Figure 14 
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APPLICATION TO PREASSEMBLED DMPC BILAYERS 


A Monte Carlo simulation using our DMPC coarse grain model is about 100 times 
more efficient than the corresponding atomistic simulation. The density profiles of the 
various groups are in very good agreement with those obtained from the corresponding 
simulation. 



Figure 15 
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SELF-ASSEMBLY OF A DMPC BILAYER 


Using the coarse grain model it takes only a few days to simulate the self-assembly of 
a DMPC bilayer from a random initial state on a typical workstation. The assembly 
process is far more efficient if we switch from a Monte Carlo simulation method to a 
molecular dynamics simulation method. 



Figure 16 
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SIMULATING LARGE SYSTEMS 


Using these models it is fairly easy to simulate large systems such as the DMPC 
system depicted below which contains 1024 phospholipids. Note that the color scheme is 
somewhat different than previously described: light blue, water; dark blue, hydrocarbons; 
yellow, hydrophilic groups. 



Figure 17 
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DMPC DIFFUSION 


These figures contain the root mean-square displacements of DMPC molecules 
during simulations employing atomistic (top) and coarse grain (bottom) models. 
Estimates for the diffusion constant, can be obtained from the limiting slopes of these 
curves. Atomistic simulation gives a value of 6.5 X 10' 8 cm 2 /sec which is in the range of 
experimental values available for this system. The coarse grain model gives a value of 
6.3 X 10' 6 cm 2 /sec or roughly 100 times larger than the atomistic simulation. This 
discrepancy is both worrying and encouraging. Worrying because it’s different than 
experiment but encouraging because it may mean that these models may evolve much 
faster than expected. We are investigating this further. 



Figure 18 
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SELF ASSEMBLY INTO A REVERSE HEXAGONAL PHASE 


As an added test we studied the self-assembly of a reverse hexagonal phase. We were 
only able to study systems large enough to contain one or two repeat units of this phase. 
However, all three studies of this phase resulted in the self assembly of the correct 
number of structures for the size of the systems studied. Below is an example from a 
system which was large enough to contain one repeat unit. 



A Phospholipid, Alkane, Water Solution 


V'i'f'f 
mmmm 
Zm mm m 
mmmm 

Simulation 



Experiment 



Figure 19 
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INITIAL STAGES OF A SIMULATION CONTAINING TWO REPEAT UNITS 


Monte Carlo simulation was used to randomly construct the initial structure and to 
relax it somewhat. However, as noted earlier for lamellar systems, Monte Carlo 
simulations were not very efficient at forming aggregates. 





Monte Carlo Simulations 


Initial 


After 8 
CPU weeks 



Figure 20 
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SELF ASSEMBLY OF A SYSTEM CONTAINING TWO REPEAT UNITS 


Switching to molecular dynamics from Monte Carlo for this simulation lead to rapid 
evolution of the system into a structure containing the two repeat units. 



Figure 21 
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A COARSE GRAIN SUCCESS STORY 


Our studies show that it is indeed possible to construct models which mimic specific 
systems fairly accurately. The range of systems covered by a consistent set of such 
models needs to be broadened. 


Course Grain Simulation Summary 

-About 10,000 times faster for phospholipids 

-Have adapted to a number of surfactant 
systems both ionic and non-ionic with both 
implicit and explicit solvation models. 

-These seem to form the basis of a self- 
consistent set of coarse grain models. 

- Huge amount of work involved in creating 
these models. 

- Still need to develop models for longer length 
and time-scales. 


Figure 22 
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MESOSCALE MODELS 


The coarse gain models that we have constructed are still too slow to study all of the 
phenomena of interest. So we must consider further simplification of the models to the 
mesoscale level. The sites are becoming quite large and there are only a few of them per 
molecule in the systems that we are studying. As a result the molecular shape, something 
which is fundamental in determining the behavior of the system that we are interested in, 
becomes distorted and constructing accurate models becomes problematic. 


Mesoscale 

DPD Models 

Can be adapted to some classes of problems 
involving specific molecules. 

e.g. polymers in a small molecule solvent. 

We encountered many cases where we could 
not adapt these models. Typically, for systems 
containing molecules of intermediate size, 
e.g. typical surfactants 


Figure 23 
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A POTENTIAL SCHEME FOR CONSTRUCTING A MESOSCALE MODEL 


Instead of using sites to represent atoms with a single molecule one could use sites to 
represent atoms across adjacent molecules as illustrated below. However, this approach 
is not entirely satisfactory in part because it restricts the association of molecules. 



Figure 24 
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SUMMARY 


Atomistic and coarse grain models can be constructed to fairly accurately mimic 
surfactant and phospholipid systems and together they provide complementary 
information. The application of mesoscale models to these systems is problematic and 
basic model construction strategies need to be explored further. 


Summary 

- Atomistic Models are fairly mature in 
many cases... still need to exercise care. 

- Have demonstrated that coarse grain 
models can, after considerable effort, be 
made to accurately describe specific systems. 

- Typical Mesoscopic models can be made to 
work in some cases. Likely need more 
complicated and accurately derived models 
to make them work in general. 


Figure 25 
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PRESENTATION OUTLINE 


A current issue facing our aging fleet of commercial and military aircraft is what is 
the residual strength of comments, such as fuselage and wing skins, after undergoing 
damage such as corrosion and fatigue crack initiation. Corrosion and fatigue crack 
initiation are inherently multiscale problems. The following discusses why we are 
looking at fatigue crack initiation and tools developed to study how fatigue cracks initiate 
and begin to propagate. These tools have been developed to study 2D representations of 
metallic polycrystals at the mesoscale that can then be used as part of a multiscale 
simulation involving continuum as well as atomics models. Recent work has involved 
transitioning these tools to 3D. 


Atif of i ah TT11 i f '1 i in 

• Overview: 

- What is the problem, and why are we working on it? 

• Mesoscale (Polycrystal) Modeling 

- Modeling Issues 

- 2D Simulation Results 

- Recent 3D Progress 

• Atomistic (Grain and Particle Boundary, Crack Tip) 
Modeling 

- Modeling Issues 

- Preliminary Simulation Results 

• Observations To Date 


Figure 2 
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MOTIVATION: FATIGUE CRACKS INITIATE IN ALUMINUM ALLOYS IN 

AGING AIRCRAFT— HOW? 


As current aircrafts are flown well beyond their original design lifetime concerns and 
research begin to focus on maintenance and repair as well as understanding how damage 
occurs. A major damage mechanism of aircraft fuselages is fatigue cracking. An aircraft 
is exposed to a corrosive environment that damages the outer surface of the aluminum 
alloy components. Fatigue cracks can then nucleate and propagate from features such as 
corrosion pits. To better see the complex physics that is going on we reduce our length 
scale to the meso-scopic level at which polycrystals become visible. Below this scale are 
networks of dislocations opening up locations for corrosive elements to penetrate deeper 
into the allow. 



Figure 3 
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MULTISCALE MODELING: THE CONCEPT 


Most analyses model components as continua. However, this approach smears out 
smaller scale features that may be the direct cause of failure at the macro-scale. Here a 
pipe is modeled with a strain field due to thermal boundary conditions on its internal 
surface. An area of high temperature can be zoomed in on and discretized to the meso- 
scopic scale, explicitly representing polycrystals. This simulation can be used to 
determine if the temperature field is resulting in decohering grain boundaries and crack 
initiation. The grain boundaries or crack tips can also be further discretized to the atomic 
scale to determine grain boundaries strength and behavior. 



Figure 4 
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POLYCRYSTAL MODELING 


Computational modeling of the mesoscopic scale at which many grains or crystals are 
explicitly represented can be broken in several steps. The steps presented here include 
determining the grain geometry and meshing for finite elements analysis, determining the 
constitutive model to be for the grains, determining the constitutive model to be used to 
describe the behavior of the grain boundaries and the grain-particle interfaces, and finally 
statistical modeling for determining grain geometry and constitutive properties. 


Poly crystal Modeling 

• Geometry and Mesh Modeling of the 
Poly crystal 

• Selecting Grain Constitutive Models 

• Selecting Grain Boundary (GB) and Particle 
Boundary Constitutive Models 

• Statistical Modeling of Geometry, 
Constitutive Properties 


Figure 5 
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MODELING: GRAIN GEOMETRY AND MATERIAL PROPERTIES 


Currently the 2D grain structure is generated using Voronoi tessellation. Each 
polygon then is taken to represent a grain. For annealed or newly crystallized metals this 
can be a good representation assuming randomly distributed initiation sites and isotropic 
grain growth at an uniform rate. This model however does not represent metals after 
processing such as rolling that results in elongated grains. The Voronoi tessellation was 
used for a first approximation. Initial simulations also did not consider sub-grain features 
such as particles, dispersoids, and precipitates. Fracture was also limited to intergranular 
debonding. Once the grain geometry is determined each grain is assigned an individual 
set of constitutive parameters from one of the models listed. Although a single 
constitutive model is chosen for the over all sample, each grain is assigned different 
parameters allowing for heterogeneity to be represented. The distribution of nano-sized 
precipitates and other hardening features are not consist from grain to grain resulting in 
heterogeneous material properties at this scale. Each parameter of the chosen constitutive 
model is randomly selected from a uniform distribution of possible values. In the case of 
orthotropic models the crystallographic orientation is sampled from the Orientation 
Distribution Function (ODF) which can be determined experimentally for the given alloy. 



* 


and Materia 





Currently modeling grain geometry using 
Voronoi tessellation 

- First excluding features such as particles, 
dispersoids, precipitates, then including 

- With intragranular cracking suppressed 

Current material models available 

- Elastic, Isotropic: Varying E 

- Elastic, Orthotropic: Varying E 1; E 2 , and b 

- Elastic-Plastic, Isotropic, Von Mises: Varying 
a yld and E 

- Elastic-Plastic, Orthotropic, Hill: Varying E,, 

a yldl’ a yld2> T yld’ aI1< ^ P 

- Properties taken front uniform distributions 

- Material angle, h, taken ifont crystallographic 
Orientation Distribution Function 



Figure 6 
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MODELING: GRAIN BOUNDARY CONSTITUTIVE RELATIONSHIP, CCZM 


The constitutive model used to describe the grain boundary response is a traction- 
displacement relationship. A coupled cohesive zone model from Tvergaard and 
Hutchinson was chosen from implementation through interface elements. Due to the 
random orientation of the grain boundaries mixed mode fracture is expected along the 
grain boundaries. For this reason a couple model was chosen. The interface elements are 
zero volume elements placed between elements in different grains. These elements were 
placed along all grain boundaries in the sample. Since the goal of the final simulations is 
to observe crack initiation the use of the cohesive zone model and interface elements 
allows for natural decohesion between grains due to the local stress and strain fields 
rather than explicit placement of a crack by the user. 


Relative Di ^placement 



fA 




MH 


Interface elements 


Modeling: Grain Boundary Constitutive 
Relationship, CCZM 

\ * Using Tvergaard-" Hutchinson 

jX Coupled Cohesive Zone Model 


tp - Strength 
k 0 ~ Initial Stiffness 


Release Rate 


• Interface Finite Elements 

- Placed between all grains 

- Follow CCZM 

No explicit introduction, of cracks. 


Interface elements allow grain boundaries to decohere to form cracks. 
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OBSERVATIONS UNDER MONOTONIC LOADING 


To observe the influences of geometry, constitutive modeling, and grain boundary 
properties on crack initiation a parametric study was conducted. The affect of geometry, 
grain constitutive model and parameters, and grain boundary parameters were studied. 
Alteration of the geometry had little affect. For all geometries generated initiation was 
primarily in Mode I. For elastic constitutive modeling, initiation paths were most 
sensitive to the presence of orthotropy. The addition of plasticity was shown to shift 
some damage to plastic deformation of the grains. For samples with large variation in the 
grain boundary strength the strain level for initiation was lowered. This was lowered 
even further when orthotrophy of the grains was introduced. The key observation was 
the role of the relationship between the grain boundary strength and the grain yield stress. 
This relationship determined how much plasticity was seen in the grains and whether or 
not decohesion occurred in the grain boundaries. 


Observations under Monotonic Loadin. 

• Geometry 

- Initiation primarily in Mode I 

• Grain Material Modeling 

- Elastic, Orthotropic 

• Initiation path sensitive to orthotropy 

- Von Mises Strain Hardening Plasticity 

• Addition of plasticity shifts some damage 
away from GB decohesion 

- Hill Perfect Plasticity 

• Absorbs most of the damage, little or no 
GB decohesion 

•GB Properties 

-Variation in GB strength reduces initiation strain level 
-When combined with orthotropic grains, reduces overall strength 
even further 

-Relationship between mean value and ranges of GB strength and grain yield 
stress most affects overall response 


Y 
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Figure 8 
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FATIGUE SIMULATION RESULT 


After completing the parametric study a fatigue test was conducted. The grain 
geometry shown was used with 100 grains modeled using Hill plasticity with mean 
primary yield stress of 500 MPa ± 5%. The mean grain boundary strength was also 500 
MPa ± 5%. The sample was then loaded under displacement control to a strain of 0.69%, 
98% of the approximate global yield strain, and unloaded twice. The third loading was 
taken to a strain of 0.70%, 99% of the approximate global yield strain. Initial decohesion 
can be seen at point 1, the first peak loading. At unloading the stress throughout the 
polycrystal did not return to 0, resulting in the presence of residual stresses at the 
beginning of the second loading. Additional opening can be seen at the second peak. 
Again residual stresses were present at the beginning of the third loading. This 
demonstrates the accumulation of damage possible throughout the sample and initiation 
due to fatigue loading. Stress concentrations can also be seen at grain boundary triple 
points at the ends of the decohering grain boundaries. At these locations the adjacent 
grain boundaries are at high angles which would require more Mode II decohesion. It is 
possible that at these points the initiating crack might propagate into the grains. 
However, the current implementation does not allow for cracks to transition from 
intergranular to intragranular. 
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Damage accumulation seen under 
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Figure 9 
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POLYCRYSTALS WITH INCLUSION PARTICLES 


After testing of the 2D polycrystal samples additional refinement was added through 
the addition of sub-grain sized particles. Again the nano-sized particles are considered to 
be smeared out and accounted for through the variation in the constitutive parameters of 
the grains. Where the nano-sized percipitates help to increase strength larger particles, or 
inclusions, can actually be detrimental to the strength of an alloy. Experimental 
observations have also shown hard particles to be initiation sites for fatigue cracks. 


Polycrystals with Inclusion Particles 

on Particles: 


Inclusion Particles: 

insoluble, undissolved or precipitated coarse 
particles 

formed and distributed heterogeneously 

Crack Initiation Observations 

• cracks form at broken inclusion particles. 

• cracks form at decohered particle-matrix 
interface of inclusions and inclusion clusters. 







..- .1" • -.:iv * 

AA 2024-T3 sheet, 500X 



Simulations 

• discretely and statistically 
represent particles 

• Perfonn simulations to 
investigate the influence of 
inclusions on fatigue crack 
initiation and propagation. 


100 jam 


Figure 10 
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MODELING: PARTICLE-MATRIX BOUNDARY CONSTITUTIVE 
REUATIONSHIP, CCZM 


Larger particles have been observed to be much harder than the surrounding matrix 
and well as brittle. These particles have also been shown to be poorly bonded to the 
matrix and thus sites for decohesion, void formation, and crack initiation. Cracking of 
the brittle particles is ignored for these simulations but decohesion is included. Again the 
particle-matrix boundary is modeled using the Tvergaard and Hutchison couple cohesive 
zone model. Interface elements are placed between the particles and the grains. Again 
this allows for natural decohesion of the particles from the matrix forming voids and 
initiating cracks. 


M 



Particle 


Matrix 

(Grain) 


Relative D isplace mein 




e™ Matrix B 

auunship, CC 

Again Using Tvergaard-Hutchinson 
Coupled Cohesive Zone Model 
t p . Strength = 0.(55 GPa 
8 C , Critical displacement - 0. 1 urn 
k 0 , Initial Stiffness = 5000 GPa/rnm 2 
G c , Critical Energy = 2.5 J/m 2 
Release Rate 


Interface Finite Elements 

- Placed between particles and matrix 

- Follow CCZM 


Interface elements 


No explicit introduction of cracks. 

Interface elements allow particle-matrix boundaries to decohere to nucleate cracks. 


Figure 1 1 
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0.4% APPLIED STRAIN 


A polycrystal sample was generated with elastic, isotropic particles with an average 
grain size of 0.017 mnr and particle volume fraction of 2%. The grains were modeled 
using Von Mises plasticity with a yield stress of 0.5 GPa. The particle-grain interfaces 
were given a strength of 0.05 GPa. The sample was then loaded under an applied strain 
of 0.4%. The sample showed stress concentrations are the particles as expected. 
Debonding of the particles was also seen. 




Figure 12 
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DEBONDED PARTICLES PROVIDE NUCLEATION SITES FOR 

MICROCRACKS 


The debonded particles provide sites for nucleation and propagation of cracks as seen 
in this example. 



Figure 13 
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STEPS FOR TRANSITIONING TO 3D 


Wlnle 2D plane strain or plane stress approximations give initial insights, the problem 
of crack propagation in metallic polycrystals is a 3D problem. The tools discussed so far 
have been for 2D representations and analyses. Transitioning theses tools to 3D will 
require several steps. First a 3D representation needs to be generated. The subsequent 
meshing of such a geometry is more difficult than for 2D samples. Again material 
properties are assigned to each grain and boundary conditions are applied. The large 
increase in the number of finite elements demands extra consideration before choosing a 
solution technique. The needed storage for such a representation and memory and 
computational time for solving will put a strain on computational resources that needs to 
be considered. 






3 



• Creating grain geometry 

• Meshing for FEM 
simulations 

• Assigning individual 
material properties for 
each grain 

• Determining appropriate 
boundary conditions 

• Selecting solution 
techniques 

• Determining necessary 
computational power 



Figure 14 
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GENERATING 3D POLYCRYSTAL MESH 


For a small example a 50 grain sample was generated using a 3D Voronoi 
tessellation. The sample was meshed using QMG 2.0 an octree mesh generator created 
by Professor Steven Vavasis at Cornell University. The resulting mesh had over 800,000 
10-noded tetrahedral elements and over 3.3 million degrees of freedom. This is a huge 
increase from the 12,000 degrees of freedom for a 100 grain 2D sample. The overall 
project of which this research is only a small part will involve several different software 
applications each operating on the same data. To facilitate the use of a single set of data 
by many application the mesh data was stored in a database using SQLServer 2000. 
Simple SQL queries can be written to extract the necessary information into a suitable 
format for the desired application. This also allows for our collaborators at other 
universities including Mississippi State University and the College of William and Mary 
to share a single set of data. 


Generating 3D Polvcrv 


50 Grain, 3D Voronoi 

Tessellation L. 

815,020 10-noded 
Tetrahedra generated by 
QMG 2.0 (Vavasis, 

Cornell) 

- Octree mesher 

- 148,605 vertices 

- 969,333 edges 

- 3,353,814 DOF 

- 280,518,175 non-zero 
matrix entries 

Data source 
( storage/ query/ retr e val ) 

SQLServer 2000 

Mesh density = number of tetrahedra sharing a given vertex 



Mesh Density 

■582 


438 


294 


150 
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TENSION TEST UNDER APPLIED DISPLACEMENT 


Each grain in the sample was randomly assigned material properties for an elastic, 
isotropic constitutive model. Three sides of the cube were then fixed as shown. The 
negative XZ face was fixed in the Y direction. The negative YZ face was fixed in the X 
direction. The XY face was fixed in the Z direction. Finally a strain of 1% was applied 
in the Z direction to the positive XY face. 


ension Test Under Applied 



• Elastic Isotropic grain 
material 

- E: randomly selected 
from uniform 
distribution with range 5 
- 1 5 * 1 0 6 

• Applied Displacement 
inZ, 1% 



Figure 16 
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TENSION TEST UNDER APPLIED DISPLACEMENT 


The problem was solved using 64 servers. Each server has 2 processors of which 
only 1 was used. Using the 64 processors the solution took less than 28 minutes. Due to 
the variation in the Young’s Modulus of each grain the resulting displacement field is not 
uniform. The wireframe polycrystal shown is the displaced shape. The Poisson effect 
can be seen as the YZ and XY faces have contracted. The color contour is shown over 
the original configuration. Visualization was created using VTK. 


Tension Test Under Applied 
Displacement 


• Solved on: 

- 64 Dell 1550 servers 

- Using 1 CPU each 

- 2x1 GHZ PHI CPU’s 

- 2GB RAM 

- Giganet interconnect 

- Windows 2000 
Advance Server 

- MPIPro 1.6.3 

• Solve time: 1664.43 s 
(< 28 min) 

- Preconditioned CG 

- Tolerance = le- 15 

- 4199 iterations 



Wireframe - displaced shape 
Displacement contour - \ \u\ \ 


Figure 17 
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COUPLING MD TO FINITE ELEMENTS 


Techniques are being investigated to transition between molecular dynamics 
representations at the atomic scale and finite element representations of the meso-scale. 
One technique currently be looked at is overlapping a region of atoms over the finite 
element representation. The atomic scale captures the nonlinear behavior in small areas 
of interest such as at a crack tip which is then added to the linear response determined by 
the finite elements. 


Coupling MD to Finite Elements 

• Instead of resolving mesh to atomic scale, 
overlap continuum region with atomistics 


Tcsi Mcsh*. quasi 2D 
Silicon cracking on I H 
piano 

^ 20 -Boded: bricks 

v— I$~j\oded wedges 
- — i/4 - point wedges 

A inm&Oe Region 
(diameTcr 7 


Layer thickness; tXR&Oi 
Width of am 


Figure 18 
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WORK IN PROGRESS 


There is still much work to be done in the areas discussed in this presentation. 
Although simple 3D poly crystal examples have been created and solved, several of the 
tools available in 2D need to be implemented in 3D. These include generating and 
meshing samples with hard particle inclusions as well as 3D cohesive elements for grain 
boundary and particle-grain interface decohesion. In order to facilitate meshing and 
solving of the samples, parallel implementations are needed. We are also still 
investigating techniques for bridging the length scales between meso- and atomic scale 
representations. These include techniques such as quasi-continuum approach and 
overlaying an atomic region as shown previously. The goal of this work is to incorporate 
the pieces shown here along with other simulators to create an adaptive software 
environment in which to study multiscale phenomenon. 


Work in Progress 

• Model polycrystal samples with particles and cohesive 
separation in 3D. 

• Bridge the length scale between meso- and atomic scales: 
quasi-continuum approach, overlaying atomic region. 

• Understand relation between atomistic modeling and CZM 
in simpler systems, e.g., single crystals. 

• Implement parallel versions of meshers and simulators: 
essential for 3D 

• Create adaptive software environments to facilitate 
multiscale simulations. 


Figure 19 
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Scale-Bridging in Dislocation-Based Multi-Scale Modeling 


Kyung-Suk Sim 
Brown University 
Providence, RI 
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INTRODUCTION 


Three examples of scale bridging problems in nano- and micro-mechanics of solids 
will be discussed in this presentation. 


Scale-bridging in dislocation-based multi-scale modeling 

Kyung-Suk Kim 
Brown University 
Providence, RI 02912 

Introduction 

1. Generalized stacking fault energy, dislocation-nucleation instability, 
and its effects on flow stress of the slip process. 

2. Multiple-dislocation cooperative slip bridging dislocation nucleation 
processes and pile-up termination processes. 

3. Continuum plastic-flow-field transition from small scale to a large 
scale, caused by elastic and plastic anisotropy incompatibilities. 

Supported by NSF/MRSEC, Ford and ONR 
Presented at NASA Langley, March 6, 2002 


Figure 1 
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CHARACTERISTICS OF DISLOCATIONS AT DIFFERENT LENGTH SCALES 


Plastic flow of crystalline solids is primarily governed by the motion of dislocations. 
The motion is controlled by various dislocation characteristics at different length scales. 
Therefore, the plastic flow and the failure processes led by plastic flow must be analyzed 
by multi-scale modeling, if responsible micro-mechanisms for an upper-scale phenomena 
are to be understood and interlinked between different length scales. This slide shows 
typical characteristics of dislocations at different length scales, that have to be modeled 
with scale bridging. 


Plasticity theory of crystalline solids 

Dislocation theory (1-10 nm, 10-100nm, 0. 1 -5 pm, 0.2-10 pm, 2-100 pm) 

Core-level dislocation model; configurational degree of freedom in the core 
Line-level dislocation model; configurational degree of freedom in the curve 
Dislocations/atmosphere interaction model; degree of freedom in the atmosphere 
Dislocation - dislocation interaction model; degree of freedom in the forest 
Dislocation - wall interaction model; degree of freedom near a wall 

Continuum theory ( 0.2-10 /uni , above 10pm) 

Strain gradient theory 

Single crystal plasticity 

Twinning and phase transformations 

Grain boundary sliding and grain growth 

Diffusion-coupled plasticity 

Void nucleation and growth 

Polycrystal plasticity and texture formation 


Figure 2 
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MAP OF DEFORMATION-MEASUREMENT TECHNIQUES 


Experimental techniques should be also scale-bridged because the resolution of an 
experimental technique has a limited domain of validity on the scale map of the field of 
view and the strain resolution. In multi-scale modeling the model at a specific length 
scale has to be checked with experimental verification at the corresponding physical 
length scale, and proper experimental techniques have to be used for the validation and 
the scale bridging. 


Map of Deformation-Measurement Techniques 

K.-S. Kim, Nano & Micromechanics Laboratory, Brown University 


Field of View (Gage Length in m) 



HRTEM High Resolution Transmission Electron Microscopy 

CFTM Computational Fourier Transform Moire 

AFM Atomic Force Microscopy 

SEM Scanning Electron Microscopy 

SRES Surface Roughness Evolution Spectroscopy 


LDLM Large Deformation Laser Moire 
FGLM Fine Grating Laser Moire 
LSI Laser Speckle Interferometry 

DIC Digital Image Correlation 


NSF Award No. CMS-0070057, Engineering Directorate (Program Manager: Dr. Jorn Larsen-Basse) 



Figure 3 
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GENERALIZED STACKING FAULT ENERGY 


Some irreversible slip processes at a nano- and micro- meter scales are predominantly 
determined by an instability of the motion of a dislocation. 

The characteristics of some slip process can be described by the generalized. 
Stacking fault energy. The generalized stacking fault energy can be obtained by ab-initio 
calculations, and it can be connected to a continuum mechanics modeling of the slip 
instability as initiation of dislocation motion. 
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SINGLE-ASPERITY CONTACT FRICTION 


The friction force required to slide a single asperity contact used to be considered to 
be proportional to the contact area, since Bowden and Tabor made the conjecture based 
on macroscopic experimental observation. The friction force divided by the contact area 
is called the contact stress, and indeed the friction stress has been measured to be constant 
within a contact radius range of 3 - 14 nm in UHV AFM friction experiments. The 
friction stress has been also observed to be constant in Surface-Force-Apparatus (SFA) 
experiments within a contact radius range of 40 - 250 microns. However, the friction 
stress of the AFM experiment is about 1/30 of the shear modulus of the contact joint, and 
that of the SFA is about 1/1300 of the shear modulus. These results contradict the 
Bowden and Tabor conjecture, and the friction stress between the two different length 
scales could not be measured so far. 



Figure 5 
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MAP OF SINGLE-ASPERITY-FRICTION MECHANISMS 


Recent analysis based on the Peierls-Rice formalism with the generalized stacking 
fault energy shows that the single asperity friction has different mechanisms of slip at 
different length scales of the contact radius. At a nano-meter scale the contact area is so 
small that it can not contain a well-defined dislocation loop and corresponding 
dislocation core structures for the slip process; therefore, the slip process is considered to 
be a concurrent slip and the friction stress is close to the theoretical strength of the joint. 
As the contact area becomes larger than a certain critical size, a mobile dislocation is 
nucleated at the edge of the contact and it self-annihilates due to its line-tension energy 
before the next dislocation is nucleated, assisting the slip process; this slip mechanism is 
called “Single Dislocation Assisted (SDA)” slip. As the contact size gets to be larger, the 
leading dislocation is stabilized with in the contact area before the trailing dislocations 
are nucleated; this mechanism is called “Multiple-Dislocation-Cooperated” slip and the 
friction stress is contact-size independent. 



J.A. Hurtado and K.- S. Kim. Proc. R. Soc. Lond. A, 455 , 3363 - 3384 (1999) 
& ibid. 3385-3400 (1999). 


Figure 6 
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CONCLUSION ON GENERALIZED STACKING FAULT ENERGY 


The evaluation of the generalized stacking fault energy is sensitive to the modeling 
scheme such as the density functional theory or an embedded-atom- potential method. It 
is a good testing entity for checking the validity of the modeling scheme. Other utilities 
of the generalized stacking fault energy as a scale-bridging scheme are listed on the slide. 


Conclusion on Generalized Stacking Fault Energy 

for degree-of-freedom reduction in a dislocation core 

Generalized Stacking Fault Energy of Aluminum (111) 



G. Lu, N. Kioussis, V. Bulatov and E. Kaxiras. Phys. Rev. B, August 2000. 

1. Generalized Stacking Fault Energy is identified as scale-bridging 
function for dislocation-core representation. 

2. Critical Experiments have been identified to understand the scale- 
dependent-mechanism transition in single-asperity-contact slip. 

3. Friction stress is found sensitive to lower-scale modeling property. 


Figure 7 
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2. MULTIPLE-DISLOCATION-COOPERATIVE SLIP MODELING (A356-T6) 


Another scale-bridging example for multi-scale modeling is discrete dislocation 
modeling for multi-dislocation cooperated slip processes. The example is the inverse size 
effect of the plastic flow stress observed in a cast aluminum A356-T6. The flow stress for 
a fine structure of secondary dendrite arm spacing is lower than that of a coarse structure. 
The cast aluminum is a multi-phase alloy containing pro-eutectic phase surrounded by 
silicon- aluminum eutectic phase. The pro-eutectic phase is softer than the eutectic phase 
as shown on the hardness distribution near the interface between the two phases. 



Figure 8 
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DISLOCATION BLOCKING MECHANISMS OF A A356-T6 


This slide shows dislocation slip bands in the pro-eutectic phase, blocked by the 
eutectic phase. This type of plastic flow requires multi-scale modeling of the collective 
dislocation motion. The first inset shows an optical view graph of the slip bands and the 
other insets show AFM view graph of the slip bands. The height of the slip steps in the 
AFM view graphs show that dislocations are nucleated in the pro-eutectic phase and 
blocked by the eutectic phase. 



Figure 9 
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DISCRETE DISLOCATION MODELING 


For the multi-scale modeling of the collective dislocation motion Needleman and his 
associates have developed a discrete-dislocation FEM methods. The method employs 
mechanism-based constitutive rules of the dislocation motion, which requires the long- 
range interaction stress field. The stress field is computed by an FEM analysis. The 
constitutive rules coupled with the FEM analysis make the multi-scale modeling 
algorithm. 


Constitutive Rules 

jj Source distribution (Strength & Position) 
Nueleation rules 

Obstacle distribution (Strength & Position) 
Mobility rules rb = Bv 

Interaction is treated by FEM. 


U » L/2T 

W as/.?. 



U = Lf 2T 



Figure 10 
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CONCLUSION ON MULTIPLE-DISLOCATION COOPERATIVE SLIP 


The conclusion page is self explanatory. The left inset shows the inverse-size effect 
observed on the stress strain relation predicted by the discrete dislocation modeling. 
Right side insets show dislocation-slip-band configurations for different-size micro- 
structures, generated by the computational modeling. In smaller microstructures the slip 
band can penetrate and glide through the eutectic phase, while the slip bands can only 
penetrate a skin depth in the eutectic phase. The skin depth forms a boundary layer of the 
dislocation penetration. 


Conclusion on Multiple-Dislocation Cooperative Slip 

with Discrete-Dislocation Modeling of Plasticity 

2 
1.5 

io a 

i 

0.5 
0 

1 . Discrete dislocation model captures the qualitative 
behavior of multi-phase plasticity at micro-scale. 

2. Inverse size effect in multi-phase has been observed 
and modeled with a discrete-dislocation model. 

3. Load re-distribution in multi -phase materials is 
found sensitive to micro-plasticity. 

Benzerga, A. A., Hong, S. S., Kim, K.-S., Needleman, A. and Van der Giessen, E., 

"Smaller is softer: an inverse size effect in a cast aluminum alloy," Acta Materia lia. Vol. 49, No. 15, pp. 3071-3083, 2001. 



Ptsuxc 3: 



□S3 


Figure 1 1 
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SYMMETRY TRANSITION FROM ASYMPTOTIC CRACK-TIP PLASTIC TO 

ELASTIC FIELDS 


The asymptotic mode-I crack-tip plastic-strain fields on (101) have been considered 
identical for [10-1] and [010] cracks in both BCC and FCC crystals, because the 
continuum yield surfaces are identical for the two crack orientations [Rice, 1987], We 
will see the experimental results on the next slide. 



Figure 12 
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EXPERIMENTAL OBSERVATIONS OF THE CRACK-TIP SLIP SECTORS 


In contrast to the predictions the experimental observations on the slide show 
completely different slip sectors for the two orientations. The upper diagrams depict 
activated slip systems predicted by Rice [1987], while the lower pictures show the slip 
traces and slip sectors in Iron-Silicon3% for the two orientations. 



Figure 13 
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COMPARISON BETWEEN EXPERIMENTS AND COMPUTATIONAL 

MODELING 


The data (a) and (b) shown in the red inset are the angular distribution of Mohr strain 
{(s xx - Syy)/2, Sxy} in the asymptotic plastic-strain fields around the crack tips of the two 
orientations. The distributions are distinctly different. However, the analytical and small- 
strain FEM simulation in (c), a large-deformation FEM simulation in (d) and discrete 
dislocation models in the other two insets could not capture the characteristics of the 
strain distribution in (b). All the computational simulations assumed isotropic elasticity. 



Figure Up Comparison bfftwwn «tpflrim«ifcaJ , FRS and hire asymptotic. anln+.inne. (a) 
AFDV'1-CUl spec.imftn at. $IJIJ ftm , Stag* 12 (b) FVS> Bp*r.inren ah Shmld and 

Kim (1991). (c) S-R asymptotic. flotation (1989) and Safl*dvaia(l991} FRS (d) C-0 
asymptotic. solution ( 1992). 


Figure 14 
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CONCLUSION ON SYMMETRIES IN CRACK TIP PLASTICITY 


Experimental measurements show that the asymptotic plastic-deformation field of (b) 
exhibits counter-intuitive unloading sectors ahead of the crack tip as shown in the first 
figure of the slide. The second figure shows the asymptotic elastic-field Mohr-strain 
distribution around the crack tips of the two orientations. The distribution shifted to the 
right is not coherent with the asymptotic plastic-deformation field, generating the 
unloading sectors. The crack-front unloading sectors make the crack growth torturous as 
shown in the third picture. The fourth figure shows the Mohr-stress distribution around 
the crack tip. As seen in this example, the scale bridging between the asymptotic plastic 
field and the outer elastic field is essential to capture proper physical processes. 


Conclusion on Symmetry Problems in Crack-Tip Plasticity 



1. Both elastic and plastic anisotropies, and associated symmetries around a 
crack tip play important roles in forming plastic-flow field. 

2. Identifying important physics is essential at any scale in multi-scale modeling 
for process-dependent characteristics such as fracture. 


Figure 15 
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CONCLUSION 


It is self-explanatory. 


Conclusion 

1 . Scale-bridging concepts such as dislocation model and cohesive-zone 
model - generalized stacking fault energy are essential for multi -scale 
modeling of plastic flow and friction. 

2. Right model has to be selected to capture right physics of the 
phenomena; discrete dislocation model captured the inverse size 
effects in multi-phase plastic flow. 

3. Identifying important physics is essential at any scale in multi-scale 
modeling for process-dependent characteristics such as fracture; 
symmetry problems are shown for single crystal crack -tip plasticity. 

4. Critical experiments are always needed at any scale for reliable multi- 
scale modeling. 


Figure 16 
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INTRODUCTION 


Yield Criteria for BCC Metals That Include 
Effects of Non-glide Stresses Based on 
Understanding the Atomic Level Features of 
Dislocation Motion 

V. Vitek 

Department of Materials Science and Engineering, 
University of Pennsylvania, Philadelphia 

J. L. Bassani 

Department of Mechanical Engineering and Applied Mechanics, 
University of Pennsylvania, Philadelphia 

SUPPORT 

NSF, Engineering Division, CMS99-00131 

ASCI through Lawrence Livennore National Laboratory 


Figure 1 
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GENERAL FRAMEWORK OF CRYSTAL PLASTICITY 


The criterion for the activation of a slip system a is that a yield function reaches a 
critical value. This function depends, in general, on the whole tensor of the applied stress 
and crystallography of the material studied. 


CRYSTAL PLASTICITY 

G ENE RAIL FRA M EWORK 

Define a yield function O for a slip system 
a which is a function of the applied stress 

a 

tensor (T 

A slip system a is activated when the yield 
function reaches a critical value 0 , 

critical 

O" [a" ; crystallography] = 0“ itical 


Figure 2 
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SCHMID’S LAW 


This criterion for the activation of a slip system states that the plastic flow 
commences when the shear stress in a given slip plane in the slip direction reaches a 
critical value. It is implicitly assumed that the components of the applied stress tensor 
other than the shear stress in the slip direction have no influence on plastic flow and that 
the sense of shearing plays no role. This criterion was deduced experimentally in studies 
of slip in hep (zinc) and fee (copper, gold) materials. 


SCHMID’S LAW 

PLASTIC FLOW WILL COMMENCE WHEN 
THE RESOLVED SHEAR STRESS ON THE 
OPERATIVE SLIP SYSTEM REACHES A 
CRITICAL VALUE 

INCLUDES TWO DISTINCT ASSERTIONS 

• The components of the applied stress 
tensor other than the shear stress in the 
slip direction have no influence on plastic 
flow 

• The critical resolved shear stress does not 
depend on the sense of shearing 


Figure 3 
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SCHMID’S LAW AND YIELD FUNCTION 


The yield function attains a very simple form when the Schmid’s law applies. For a 
slip system a, corresponding to a slip plane and a slip direction, it is just equal to the 
shear stress in the slip plane in the slip direction, i.e. to the Schmid stress. The yield 
criterion is then that the Schmid stress reaches a critical value. 


SCHMID BEHAVIOR 

The yield function for the slip system a is equal to 
the Schmid stress for this system, T 

When the critical resolved shear stress for the slip 

a 

system a is ^critical the yield criterion is 





a 

critical 


Figure 4 
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SLIP CHARACTERISTICS OF FCC MATERIALS 


The most important characteristics of the slip in pure fee metals, and also pure hep 
metals when the slip is along basal planes, are very weak temperature dependence of the 
flow and yield stress indicating no significant lattice friction Peierls stress. Close packed 
planes, i.e. {111} in fee and basal in hep, are the only slip planes and no asymmetries 
between tension and compression are observed. The Schmid’s law applies for this type 
of gliding and the only important stress component governing the plastic flow in fee 
materials is the shear stress in {111} planes in <11 0> directions. However, deviations 
from the Schmid’s law occur in the case of cross slip. 


SCHMID’S LAW APPLIES IN CLOSE PACKED 
STRUCTURES SUCH AS FCC AND HCP 

CHARACTERISTICS OF SLIP IN FCC MATERIALS 

• Well defined slip planes {111} 

• Weak dependence of the yield and flow stress on 
temperature 

• No significant lattice friction, i. e. Peierls stress 

• The same behavior in tension and compression 


Figure 5 
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SLIP CHARACTERISTICS OF BCC MATERIALS 


The most important characteristics of the slip in pure bcc metals are the very strong 
temperature dependence of the flow and yield stress indicating large lattice friction 
Peierls stress, complex slip geometry, pronounced asymmetry between tension and 
compression. The Schmid’s law does not apply. Fore reviews see [1-3], 

1. Christian, J. W ,,Metall. Trans. A Vol. 14, 1983, p.1237. 

2. Vitek, V., Dislocations and Properties of Real Materials, edited by Lorretto, M., 
(London: The Institute of Metals) 1985, p. 30. 

3. Duesbery, M. S., Dislocations in Solids, edited by Nabarro, F. R. N., (Amsterdam: 
North Holland) Vol. 8, 1989, p. 67 


CHARACTERISTICS OF SLIP IN 
BODY-CENTERED CUBIC METALS 

• Strong temperature and strain rate dependence 
of the flow stress and slip geometry 

As T •“+ OK the flow stress ~10‘ 2 G 

• Strong and unusual orientation dependencies 
and asymmetries of the flow stress and slip 
geometry 

• Slip planes are not always well defined 

• SCHMID’S LAW IS NOT VALID 


Figure 6 
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OTHER MATERIALS WITH ANALOGOUS SLIP PROPERTIES 


Materials with bcc structure are no exception. The Schmid law is not valid and strong 
temperature and orientation dependencies of the yield stress are observed in many 
intermetallic compounds and non-cubic structures. Rather, fee materials and hep gliding 
on basal planes are a special case when Schmid’s law applies. An example is the 
molybdenum di-silicide (MoS^). In this material dislocation glide occurs on are several 
slip systems. The most important is {013} <3 31] which displays large asymmetry, in 
particular for the orientation of the compression axis close to [001], 


MAIN OPERATING SLIP 
SYSTEM 

{013)<331J 

YIELD STRESS IS STRONGLY 
DEPENDENT ON ORIENTATION 

For [001] orientation of the 
compression axis MoSii is brittle 

(013)<331] has Schmid factor 
relatively large but does not operate 

But, it starts to operate for small 
deviations of the compression axis 
away from [001] 


Figure 7 
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SIGNIFICANT CONSEQUENCES OF THE BREAK DOWN OF THE SCHMID’S 
LAW FOR THE PLASTIC BEHAVIOR 


When the Schmid law does not apply the plastic deformation displays asymmetries 
that are most frequently characterized by different behavior in tension and compression. 
Such unique stress-state dependence of flow properties arising when Schmid’s law does 
not apply can significantly affect critical phenomena such as strain localization, evolution 
of the texture of grains, and failure mechanisms at mesoscopic and macroscopic length 
scales. 


CONSEQUENCE OF NON- 
SCHMID BEHAVIOR 

• Tension-compression asymmetry in 
single and polycrystals 

• Elevation of stresses due to 
restriction of slip in polycrystals 

• Strain localization - formation of shear 
bands 

• Affects evolution of textures and 
failure mechanisms 


Figure 8 
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NON-PLANAR DISLOCATION CORES 


The most common reason for the very high Peierls stress and associated asymmetries 
of dislocation motion that lead to the break down of the Schmid law are dislocation cores 
spread into several non-parallel planes. Schematic picture of such dislocation cores is 
shown below. In bcc metals the screw dislocations possess such cores. 
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SCREW DISLOCATIONS IN BCC METALS 


It is now generally accepted that in bcc metals screw dislocations control their 
characteristic plastic properties, in particular temperature and orientation dependencies of 
the yield stress as well as their slip geometry [1-3], 

1. Christian, J. W ,,Metall. Trans. A Vol. 14, 1983, p.1237. 

2. Vitek, V., Dislocations and Properties of Real Materials, edited by Lorretto, M., 
(London: The Institute of Metals) 1985, p. 30. 

3. Duesbery, M. S., Dislocations in Solids, edited by Nabarro, F. R. N., (Amsterdam: 
North Holland) Vol. 8, 1989, p. 67 


1/2[111] SCREW DISLOCATION 

• NON PLANAR CORE SPREADS INTO 
SEVERAL PLANES OF [111] ZONE 

• INTRINSIC ASYMMETRIES 

• COMPLEX RESPONSE TO APPLIED 
STRESSES 


Figure 10 
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ATOMISTIC COMPUTER MODELING OF SCREW DISLOCATION IN 

TRANSITION BCC METALS 


Computer simulation of the dislocation core was done using a molecular statics 
method and central-force many-body potentials for Mo and Ta. The core structure is 
represented using the method of differential displacements. The atomic arrangement is 
shown in the projection perpendicular to the direction of the dislocation line ([111]) and 
circles represent atoms within one period without distinguishing their positions in three 
successive (111) planes. The relative displacement of the neighboring atoms produced 
by the dislocation in the [111] direction is depicted as an arrow between them. In the 
case of Mo the core spreads in a three-fold manner into the three {110} planes while in 
the case of Ta it spreads in a six-fold manner into the same planes [1,2], 

1. Duesbery, M. S. and Vitek, V., Acta Mater. Vol. 46, 1998, p. 1481. 

2. Ito, K. and Vitek, V., Philos. Mag. Vol. 81, 2001, p. 1387. 


CORE STRUCTURE OF THE 1/2|1 1 1] 
SCREW DISLOCATION 



Core structure 
in molybdenum 



ORIENTATION 
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Core structure 
in tantalum 


Figure 1 1 
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TRANSFORMATION OF THE CORE PRIOR TO THE DISLOCATION 

MOTION 


The most important aspect of non-planar dislocation cores is that their structure 
transforms under the influence of applied stresses prior to the dislocation motion. An 
example presented here shows changes in the core of the 1/2[111] screw dislocation 
under the effect of pure shear in the [1 1 1] direction in the plane. 
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DEPENDENCE OF THE CRITICAL RESOLVED SHEAR STRESS (CRSS) ON 
THE ORIENTATION OF THE PLANE OF THE MAXIMUM RESOLVED 

SHEAR STRESS (MRSSP) 


The calculated dependence of the CRSS on the orientation of the MRSSP in 
molybdenum shows the break down of the Schmid’s law and asymmetric behavior of this 
dependence. 



Figure 13 
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TENSILE AND COMPRESSIVE LOADING 


Atomistic modeling of the dislocation motion was performed for a number of 
different tensile and compressive axes marked in the standard triangle. The three groups 
of the tensile and compressive axes correspond to three different maximum resolved 
shear stress planes. 


( Tension/Compression axesT) 


Compression 



Figure 14 
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MOLYBDENUM TENSION AND COMPRESSION 



Figure 15 
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EFFECT OF THE SHEAR STRESS PERPENDICULAR TO THE BURGERS 

VECTOR ON CRSS 


The stress component affecting the glide in the case of tensile/compressive loading 
was identified as the shear stress perpendicular to the Burgers vector, i.e. the stress that 
does not exert the Peach-Koehler force on the dislocation. Computer modeling 
combining this shear stress with the glide shear stress parallel to the Burgers vector 
demonstrates the dependence of the CRSS on this non- glide shear stress and shows that 
tensile/compressive data fit into this dependence. 
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EFFECT OF THE SHEAR STRESS PERPENDICULAR TO THE BURGERS 
VECTOR ON THE DISLOCATION CORE 


The reason why the non-glide shear stress perpendicular to the Burgers vector affects 
significantly the CRSS is that it induces changes in the dislocation core that may either 
help or hinder the eventual glide, depending on the orientation of the non-glide stress. 


EFFECT OF THE SHEAR STRESS 
PERPENDICULAR TO THE BURGERS ON THE 
CORE OF THE SCREW DISLOCATION IN 
MOLYBDENUM 



x = -0.052C 44 x = 0.052C 44 


Figure 17 
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MOST IMPORTANT PHYSICAL FINDING BASED ON THE ATOMISTIC 
MODELING OF THE CORE OF THE SCREW DISLOCATION 


The glide of the dislocation is affected and governed by shear stresses both parallel 
and perpendicular to the Burgers vector that act in the {110} planes of the [111] zone. 


PHYSICAL CONCLUSION 


THE YIELD FUNCTION DEPENDS ON 

• SHEAR STRESSES PARALLEL TO THE BURGERS 
VECTOR ACTING IN THREE {110} PLANES OF 
THE [111] ZONE 

• SHEAR STRESSES PERPENDICULAR TO THE 
BURGERS VECTOR ACTING IN THREE {110} 
PLANES OF THE [111] ZONE 

• BOTH THESE DEPENDENCIES HAVE BEEN 
DETERMINED BY ATOMISTIC MODELING OF 
THE DISLOCATION MOTION 


Figure 18 


260 






COMBINING THE RESULTS OF ATOMISTIC STUDIES WITH CONTINUUM 

CRYSTAL PLASTICITY 


Based on the findings of atomistic modeling of the glide of 1/2<1 1 1> dislocations in 
Mo and Ta the yield function depends on shear stresses in the {110} planes both parallel 
and perpendicular to the Burgers vector. A feasible approximation is that the yield 
function is a linear combination of the Schmid stress and non-glide stresses. The 
corresponding coefficients need to be determined so as to reproduce the results of 
atomistic modeling [1-3], 

1. Qin, Q. and Bassani, J. L., J. Mech. Phys. Sol. Vol. 40, 1992, p. 813, 835. 

2. Bassani, J. h.,Adv. Appl. Mech. Vol. 30, 1994, p. 191. 

3. Bassani, J. L, Ito, K. and Vitek, V., Mat. Set. Eng. A Vol. 319, 2001, p. 97. 


YIELD FUNCTIONS WITH NON-GLIDE STRESSES 

stress tensor 


- , * , , a a ^ a 

timid stress on slip system a: x = n • a • m 


slip plane normal 


slip direction 



components: 




• a* m 
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”nand m“ are crystallographic vectors that resolve 
each of the ( rj =1, n /; „) non-glide stress components 
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non-glide stress yield parameters 


Figure 19 
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YIELD FUNCTION FOR APPLICATION OF PURE SHEAR STRESSES 
PARALLEL TO THE BURGERS VECTOR 


When a shear stress parallel to the Burgers vector is applied in a MRSSP making 
angle % with the slip plane we can choose the shear stress in any other plane of the [111] 
zone as the non-glide stress entering the yield function. 
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YIELD FUNCTION FOR APPLICATION OF THE PURE SHEAR STRESS 
PARALLEL TO THE BURGERS VECTOR 


The atomistic studies suggest to choose the non-glide shear stress entering the yield 
function the shear stress in one of the two {110} planes that are not the slip planes. Here 
we choose the shear stress in the plane. 


Choosing the non-glide stress as the shear stress in the (110) plane 
T ng l0) = X MRSS cos 0c + 60) and the yield criterion reads 

"^Schmid + mV = T v,RSS [ C0S < 7.) + «COS(X + 60)]= T cr 

CRSS = x cr [cos(x) + «cos(x + 60)] 



X(DEGREES) 

Figure 21 
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GLIDE ON ALL POSSIBLE SLIP SYSTEM IN A BCC CRYSTAL 


There are twenty four <1 1 1>{ 1 10} slip systems (considering positive and negative 
slip directions as different). Each of these slip systems follows the same yield criterion 
but they all may operate at the same time. 
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TAYLOR MODEL WITH NON-GLIDE STRESSES 


For the polycrystal made of randomly oriented grains of bcc Mo we assume that the 
strain in each crystal is the same as the macroscopic strain (Taylor hypothesis). At the 
same time each grain satisfies the yield criterion involving the non-glide stresses. The 
Taylor factor demonstrates the tension-compression asymmetry in this poly crystal 
induced by the effect of non-glide stresses. 


Taylor Modal with Hon-Glida Stresses 

Consider a polycrystal of randomly oriented BCC grains each satisfying the yield 
criteria: z* a = r a +arf = r*“ ■ Neglecting elastic strains and assuming that the strain in 
each crystal is the same as the macroscopic strain (Taylor hypothesis), a quadratic 
programming problem is used to solve for the minimum of 5 slips in each crystal, which 
gives an upper bound to the limit yield surface. For Schmid behavior («=0) the classical 
Taylor factor is 3.07 times the slip-system y ield stress in tension and compression. 



& C 

&t+V C 



a a 

Taylor factor for a random polycrystal 
with T* a = T a +arf at the grain level 


Figure 23 
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EFFECT OF NON-GLIDE STRESSES ON THE YIELD SURFACE AND 

HARDENING 


The calculated yield surface and stress-strain curves demonstrate asymmetries 
induced by the effect of non-glide stresses and thus the effect of non-Schmid behavior 
that occurs on the macroscopic scale owing due to the atomic level properties of 
dislocations in bcc metals. 



Plane stress yield surfaces for a random 
polycrystal 

Inner locus: bcc model with <2-0.6 the 

Outer locus: Schmid behavior with a= 0 

Stresses are normalized by the critical 
stress on each slip system, t 0 . 






Figure 24 
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CONCLUSIONS 


Atomistic modeling of dislocation cores and their response to applied stresses 
determine for a given slip system: 

(i) Those components of the stress tensor that drive and affect the dislocation motion. 

(ii) Relative importance of different stress components in determination of the yield 
stress. 

The yield function for crystal plasticity is defined based on these results of atomistic 
calculations. 

This yield function is then employed in studies of plastic yielding in both single and 
polycrystals. 

Atomic level behavior of dislocations percolates into the mesoscopic and 
macroscopic plastic behavior. 


CONCLUSIONS 

Atomistic modeling of dislocation cores and their response 
to applied stresses determine for a given slip system: 

(i) Those components of the stress tensor that drive and 
affect the dislocation motion 

(ii) Relative importance of different stress components in 
determination of the yield stress 

The yield function for crystal plasticity is defined based on 
these results of atomistic calculations 

This yield function is then employed in studies of plastic 
yielding in both single and polycrystals 

ATOMIC LEVEL BEHAVIOR OF DISLOCATIONS 
PERCOLATES INTO THE MESOSCOPIC AND 
MACROSCOPIC PLASTIC BEHAVIOR 


Figure 25 
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FUTURE DEVELOPMENT 


Quantum mechanics based description of atomic interactions in transition metals and 
alloys reflecting mixed metallic and covalent bonding. 

Present development: Bond-order potentials. 

Temperature and strain rate dependence of the yield stress based on mesoscopic 
models that incorporate the results of atomistic studies via the dependence of the 
activation enthalpy on the relevant components of the stress tensor. 

Incorporation of temperature and strain rate effects together with hardening into yield 
criteria that will be used continuum studies, for example finite elements calculations, of 
deformation under complex applied stresses. 


FUTURE DEVELOPMENT 

Quantum mechanics based description of atomic interactions 
in transition metals and alloys reflecting mixed metallic and 
covalent bonding 

Present development: Bond-order potentials 

Temperature and strain rate dependence of the yield stress 
based on mesoscopic models that incorporate the results of 
atomistic studies via the dependence of the activation enthalpy 
on the relevant components of the stress tensor 

Incorporation of temperature and strain rate effects together 
with hardening into yield criteria that will be used continuum 
studies, for example finite elements calculations, of 
deformation under complex applied stresses 


Figure 26 
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PLASTICITY 

From Nano-scale to Macro-scale 
through ad hoc Homogenization 

Constitutive Constitutive Constitutive 

Model Model Model 

Lattice Scale (Nano) Crystal (Micro) Polycrystal 
(Meso/Macro) 


CEAM/UCS 


Figure 3 


Approaches 

• Continuum 

• Direct numerical simulations with various 
levels of approximations: 

-Atomistic => Macro 

• Mechanism modeling at several length scales 
with transitional (ad hoc) homogenization 

- Dislocations 

- Slip 

- Crystal 

- Polycrystal 

CEAM/UCSC 


Figure 4 
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Length-scale as Internal-state Variable 

* Strain gradient exists at all scales: 

- around dislocation lines 

- near slip planes 

- within each crystal 

- over many crystals 

- at continuum scale around geometric & material discontinuities 

* Strain fields are obtained by solving initial-boundary- 
value problems at relevant length scale 

* Constitutive relations must reflect essential structural 
characteristics at relevant length scales 

* Dislocation density provides a natural length scale in 
continuum plasticity 

- (average dislocation density) ' 12 = length scale 

CEAM/UCS 

Figure 6 
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Figure 8 
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Figure 1 1 



Figure 12 
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WHAT IS NEEDED IN THE SCIENCE OF 
DEFORMATION MODELING 

• Experimental Techniques & Tools 

• Nano/Microstructure-Properties 
Relations: Mechanisms 
Identification 

• Physics-based Nano- Micro- & 
Macro-mechanical Modeling 

CEAM/UCS 

Figure 13 

Experimental Tools 

• Recovery Hopkinson Technique 

* Compression Tests 

* Tension Tests 

* Compression-Tension Tests 

* Tension-Compression Tests 

• Isothermal and Adiabatic Tests 

* Strain Rate to 50,000/s 

* Temperature Range -200 to 1500 C 

• Tri-axial Hopkinson Techniques 

CEAM/UCS 

Figure 14 
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PLASTICITY 

Recovery Hopkinson Experiments 

Temperature: Liquid Nitrogen to 1000°C 



Figure 16 
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Equivalent Stress (MPa) 



Flow Stress 

Dislocation Motion Leads to Plastic Flow 

• Resistance to Dislocation Motion Leads to 
Flow Stress and Workhardening 

Barriers to Dislocation Motion 

• Short-range Barriers 

• Long-range Barriers 

• Viscous Drag 

CEAM/UCS 

Figure 18 
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Barriers to Dislocation Motion 


Short-range Barriers 

• Lattice Resistance = Peierls Barrier 

• Dislocations Intersecting Slip Planes 

• Solute Atoms 

• Vacancies and Interstitials 

• Precipitates, etc. 

Long-range Barriers 

• Grain Boundaries 

• Farfield Forests of Dislocations 

• Other Farfield Inhomogeneities 


CEAM/UCS 




Figure 20 
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PLASTICITY 
BCC Crystals 


Short-range Barrier 


Peierls Barrier 


Double Kink 
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Figure 21 


PLASTICITY 
Flow Stress 

T = T* + T a +T d 

x* = thermally activated part: short-range 
x a + x d = athermal part (long-range) + viscous part 

Peierls & Other Barriers 

AF = AG-x*V* 

V* = b XI, 1 = lattice spacing (bcc) 

/ = disloc. spacing (fee) 

o) = co 0 exp { -AG/kT} 

® 0 = 10 U /s to 10 13 /s 

CEAM/UCSd 


Figure 22 
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Direct Experimental Strategy 


• At fixed strain rate, obtain high -temperature 
strain-stress reiation => athermal part 

• At fixed strain, obtain high-temperature 
stress/strain-rate relation => drag part 

• Subtract athermal & drag parts from total 
flow stress => thermally activated part 

CEAM/UCSC 

Figure 23 
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True Stress,? {MPa} 





0.00 0.10 0.20 0.30 0.40 0.50 


Flow stress for indicated strain rate 
at an initial temperature of 296K 


Effect of temperature jump from 77 to 400K, 
on flow stress of niobium at 10 3 s 1 strain rate 




Effect of strain rate jump from 10~ 3 to 3,300 s' 1 
on flow stress of niobium at initial temperature of 296K 


Limiting values of flow stress at large 
temperatures, as a function of strain 


CEAM/UCST 
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'rue Stress (MPa) 



Figure 28 
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True Stress (MPa) 



T rue Strain 


CEAM/UCSm 


Figure 29 



T rue Strain 


Comparison of model predictions with experimental results for 
indicated initial temperatures and indicated strain rates 

CEAM/UCSC 


Figure 30 
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True Stress (MPa) True Stress (MPa) 



CP MO: Stress-strain Curves 


CEAM/UCSd 


Figure 31 
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CP MO: Stress-strain Curves 


CEAM/UCSq 


Figure 32 



286 





True Stress (MPa) True Stress (MPa) 



True Strain 


CP MO: Stress-strain Curves 


CEAM/UCSC 


Figure 33 



T rue Strain 


CP MO: Stress-strain Curves 


CEAM/UCSC 


Figure 34 
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True Stress (MPa) I True Stress (MPa) 



T ru e Strai n 


CP MO: Stress-strain Curves 


CEAM/UCSC 


Figure 35 



CP MO: Stress-strain Curves ceam/ucsJ 


Figure 36 
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True Stress (MPa) | | True Stress (MPa) 




Figure 38 
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True Stress (MPa) 




0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 

True Strain 

Comparison of model predictions with experimental results at a 
strain rate of 3,500/s CEAM/UCSI 

Figure 40 
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True Stress (MPa) 




0 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 


0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 


Fig. 24. Comparison of model predictions with 
experimental results at a strain rate of 3,500/s 


Fig. 25. Comparison of model predictions with 
experimental results at a strain rate of 3,500/s 



AL6-XN.0 001/s 
. Dashed Curves: Experiment 
; Solid Curves: Model Predict 


0.10 0.20 0.30 0.40 0.50 0.60 


Fig. 26. Comparison of model predictions with 
experimental results at a strain rate of 8,000/s 


Fig. 27. Comparison of model predictions with 
experimental results at a strain rate of 0.001/s 


CEAM/UCS 
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y - y r exp( -AG / kT ) 
Y r ~ b P m ® 0 ^ 
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Table 1 Constitutive Parameters for indicated Commercially Pure Meta! 



Table 2. Constitutive Parameters for Indicated Metal 



*NA= Not available 


CEAM/UCS 


Figure 47 



Constitutive Parameters for 
Various Titanium 6AI-4V 

V alues of various constitutive parameters for the three Ti64 materials 


Parameter 

P 

q 

k/G 0 (xlO _5 K _ 

‘) toCxlO 10 ^ 1 ) 

ao 


n 

Commercial Ti64 

1 

2 

6.2 

1 32 

2.4 

1560 685 

0J05 

RS-Mil-Hip Ti64 

1 

2 

6.2 

1.32 

2.4 

1900 710 

0.03 

RS-Hip Ti.S4 

1 

2 

6.2 

1.32 

2.4 

1620 680 

004 
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Single Crystal Deformation 

Crystallographic slip leaves the lattice 
unaffected 

Plastic flow consists of simple shearing due 
to dislocation motion associated with a finite 
number of active slip systems, resulting in 
both plastic deformation and rigid-body 
rotation of the material relative to the lattice 

Elastic deformation is by lattice distoii^^ M/UCS[: 


Figure 49 
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A Simple Continuum Model 

r = L: ( D-f^ 

L = Elasticity Tensor 
D = Deformatio n Rate Tensor 

CEAM/UCSC 

Figure 51 


General Continuum Model 

o .80 

t = L : (D - y ——) 

or 

0 = Flow Potential 


CEAM/UCSq 


Figure 52 
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Model vs Exprt: Annealed OFHC Copper, 4000/s 



Equivalent Strain (%) Equivalent Strain (%) 


CEAM/UCS 
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INTRODUCTION 
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Nam^rysiaMne TMo. Fliara* 


Dorel Moldovaa Dieter Wolf, Simon R. Phillpot 
and Andrew J. Ila si am 

Materials Science Division, Argonne National Laboratory 
Aigonne, IL 60439 


*work supported by the U. S. Department, of Energy 





Figure 1 
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ABSTRACT 


We have combined molecular-dynamics simulations with mesoscopic simulations to 
study the mechanism and kinetics of grain growth in a thin film of nanocrystalline 
palladium with a columnar grain structure. The conventional picture is that grain growth 
results solely from the migration of grain boundaries in response to the driving force 
associated with the reduction in the grain-boundary area. However, our molecular- 
dynamics simulations suggest that, at least in a nanocrystalline microstructure, grain 
rotations play an equally important role. Based on this insight we have developed a 
kinetic model for grain-boundary-diffusion accommodated grain rotation, extending the 
formalism of Raj and Ashby (1971) for grain-boundary sliding. We have incorporated 
this model into a mesoscopic-simulation code based on the Needleman-Rice (1980) 
variational formalism for the dissipated power. The presence of both grain-boundary 
migration and grain rotation introduces a physical length scale, Rc, into the system. The 
growth process is characterized by two regimes: if the average grain size is smaller than 
Rc, grain growth is grain-rotation dominated; by contrast, if the average grain size is 
greater than Rc, grain growth is dominated by grain-boundary migration. Our study 
reveals that the growth exponents characterizing the power-law time dependence of the 
average grain size are different for the two growth regimes. We discuss how this 
methodology can be extended to allow the study of deformation process. This 
combination of atomic-level with mesoscopic simulations then enables the investigation, 
in a physically realistic manner, of grain growth in systems containing a large number of 
grains and over long times. 


*Work supported by the U.S. Department of Energy, Office of Science under Contract 
W-31-109-Eng-38. 


302 



OUTLINE 


Atomistic simulations 

• Evidence of grain-rotation during grain growth 

Mesoscopic simulations 

• The model 

• Variational formulation for dissipated power 

• Atomistic and mesoscale modeling of a 25 grain system 

• Validation of the mesoscale model 

• Isotropic and anisotropic grain growth 

• Grain growth by grain rotation only 

• Grain growth in the presence of both GB migration and grain rotation 
Conclusions 


qjSIM: 

• Atomistic simulations 

• Evidence of grain-rotation during grain growth 

• Mesoscopic simulations 

• The model 

• Variational formulation for dissipated power 

• Atomistic and mesoscale modeling of a 25 grain system 

• Validation of the mesoscale model 

• Isotropic and anisotropic grain growth 

• Grain growth by grain rotation only 

• Grain growth in the presence of both GB migration and grain rotation 

• Conclusions 
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DETAILS ON THE MOLECULAR DYNAMICS SIMULATIONS 


An embedded atom method potential parameterized for Pd is used in these molecular 
dynamics (MD) studies. A textured columnar microstructure was selected with the <001> 
direction as the common texture axis; therefore all GBs are <001> tilt boundaries. The 
initial microstructure contains 25 grains (388,656 atoms) with an average grain size of 
about 14nm. The constant pressure MD simulations with 3D-periodic border conditions 
were carried out on this microstructure using 5th order Gear predictor-corrector algorithm 
to solve the equations of motion. The simulation temperature was about 0.9 of the 
melting temperature for the simulated potential and the simulation was carried on for 
about 7.2ns (1.4 million MD steps). 

(A.J. Haslam, et al.. Materials Science and Engineering A3 18,293-3 12, 2001 ) 


Details cm the Molecular Dynamics Simulations: 

An embedded atom method potential parameterized for Pd is used in these 
molecular dynamics (MD) studies. A textured columnar microstructure was 
selected with the <00 1> direction as the common texture axis; therefore all GBs 
are <00 1> tilt boundaries. The initial micro structure contains 25 grains (388,656 
atoms) with an average grain size of about 14nm. The constant pressure MD 
simulations with 3D-periodic border conditions were carried out on this 
microstructure using 5th order Gear predictor-corrector algorithm to solve the 
equations of motion. The simulation temperature was about 0.9 of the melting 
temperature for the simulated potential and the simulation was carried on for 
about 7.2ns (1.4 million MD steps). 

(A.J. Haslam, et al, Materials Science and Engineering A3 18,293-3 12, 2001 ) 



Figure 3 
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GRAIN GROWTH ON MD TIME SCALE 


The positions of the miscoordinated atoms in one of the (001) planes of the 
simulation cell at the beginning and end of the simulation are shown. It is obvious the 
substantial amount of grain growth that has taken place. 


Grain growth on MD time scale 

The positions of the miscoordinated atoms in one of the (001) planes of the 
simulation cell at the beginning and end of the simulation are shown. It is obvious 
the substantial amount of grain growth that has taken place. 



Aiygcmne National Laboratory 


Figure 4 
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MD SIMULATIONS 


Our MD simulations shows that there are two growth mechanisms which lead to the 
final microstructure. These are grain boundary migration and grain rotation - grain 
coalescence. The two figures bellow focusing on three representative grains in the 
microstructure by plotting all the perfect fee coordinated atoms shows a grain rotation 
grain coalescence event. 


Our MD simulations shows that there are two growth mechanisms which lead 
to the final microstructure. These are grain boundary migration and grain 
rotation - grain coalescence. The two figures bellow focusing on three 
representative grains in the microstructure by plotting all the perfect fee 
coordinated atoms shows a grain rotation grain coalescence event. 



Figure 5 
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SUMMARY OF MD SIMUALTIONS 


We have observed two mechanisms for grain growth: grain rotation / grain 
coalescence grain-boundary (GB) migration 

At the nanoscale, grain rotation and GB migration are intricately coupled and 

play equally important roles in the grain-growth process. 

MD provides input for mesoscale simulations 
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MESOSCOPIC MODEL OF A 2D POLYCRYSTALLINE MICROSTRUCTURE 


At the mesoscale level the microstructure consists of interconnected polygons 
initially generated by a Voronoi construction with periodic border conditions. In 
our representation each grain boundary connecting two triple-points is 
represented by a set of straight segments. 

GB energy and mobility 



Figure 7 
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GRAIN-BOUNDARY ENERGY 


Grain-boundary energy parameters used in the mesoscale simulations are fitted from 
MD simulations (Pd , <00 1> tilt boundaries). 


Grain-boundary energy parameters used in the mesoscale simulations are fitted 
from MD simulations (Pd , <00 1> tilt boundaries) 



(Extended Read-Shockley) 

■y(9) = sin(26>)^ - ln(sin(20))J 

— - = 1.01 Jnf 1 

b 

^•=0.70 Jm 2 

b 

(A.J. Haslam, et al.. Mater. Sd. and Eng. A3 16, 293, 2001 ) 
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Figure 8 
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THEORY OF DIFFUSION ACCOMMODATED GRAIN-ROTATION 


(D. Moldovan, et cil ., Acta Mater. 49, 3521, 2001; based on Raj&Ashby 1971) 

Based on the theory of diffusion-accommodated grain-boundary sliding by Raj and 
Ashby we have developed a kinetic theory of grain rotation in columnar polycrystalline 
microstructures. The driving force for rotation of any given grain is the aggregate torque 
on the grain due to the misorientation dependence of the grain boundaries delimiting the 
grain. 


Theory of diffusion accommodated grain-rotation 

(D. Moldovan, et tr/.. Act a Mater. 49, 3521, 2001; based on Raj&Ashby 1971) 

Based on the theory of diffusion-accommodated grain-boundary sliding by Raj 
and Ashby we have developed a kinetic theory of grain rotation in columnar 
polycrystalline micro structures. The driving force for rotation of any given grain is 
the aggregate torque on the grain due to the misorientation dependence of the 
grain boundaries delimiting the grain. 



CO — M(R) T viscous-like grain rotation 


(M(R) -rotational mobility) 


C(Q,SD GE ,T,shape) 
M(K) - p 



for GB diffusion 
for volume diffusion 


For a regular hexagonal grain: M( R) 
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Figure 9 
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VARIATIONAL FORMULATION 


As the formal basis for our mesoscale simulations we adopt the theoretical approach 
of Needleman and Rice based on variational principle for dissipative systems, originally 
formulated for GB and surface diffusion in the context of GB void growth studies. Later 
this was adapted for the simulation of curvature-driven grain growth by Cocks and Gill; 
their modification describes the rate of power dissipation due to the competition between 
the reduction of GB energy and the viscous drag during GB migration. 

(Needleman & Rice, Acta Metall., 28, 1315, 1980); Cocks& Gill, Acta Mater., 44, 
4765, 1996) 

Our extension includes an additional dissipative process due to grain rotation. 
(Moldovan, D., Wolf D., Phillpot S.R. and Haslam A.J., Acta Mater, (in press) 2002) 


Yari ati onal Formulation 

As the formal basis for our mesoscale simulations we adopt the theoretical approach of 
Needleman and Rice based on variational principle for dissipative systems, originally 
formulated for GB and surface diffusion in the context of GB void growth studies. Later this 
was adapted for the simulation of curvature-driven grain growth by Cocks and Gill; their 
modification describes the rate of power dissipation due to the competition between the 
reduction of GB energy and the viscous drag during GB migration. 

(Needleman&Rice,>fctaMeta//.,28, 1315, 1980); Cocks& Gill, Acta Mater.,44, 4765, 1996) 
Our extension includes an additional dissipative process due to grain rotation. 

(Moldovan, D., Wolf D., Phillpot S.R. and Haslam A.J., Acta Mater, (in press) 2002) 

Two dissipative processes are considered: 

• grain-boundary migration ({v},{r}) 

• grain rotation({ro },{()>}) 

n(M,M;M, {<!>}) 

The problem is reduced to finding the variational 
parameters! v'!and{o>} such that: 

5fl=0 





Figure 10 


311 






VARIATIONAL FUNCTIONAL TERM FOR A GB SEGMENT 


Here is the explicit relation of the variational functional term corresponding to a GB 
segment of length L; moving with the velocities vu and v^. 



Figure 1 1 
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VARIATIONAL FUNCTIONAL TERM FOR A ROTATING GRAIN 


The variational functional term corresponding to a rotating grain with the angular 
velocity co . The rotational mobility of the grain is M. 



Figure 12 
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ALL TOGETHER: THE GENERAL FORM 
OF THE VARIATIONAL FUNCTIONAL 


The general form of the variational functional for a microstructure comprised of 
NGgsegm GB segments and N gr ams grains. 


All Toge ther: The General Form of t he Variational Functional 

The general form of the variational functional for a micro structure comprised of 
N OBseg m GB segments and N grajns grains. 


QBsegm j N grains ,, 2 

£ y t (y n •*„ + v l 2 -s l2 )+ ^-[(v") 2 + (v" 2 ) 2 + (v"v" )] + Z Tj<Dj+ ~2M 


Find the variational parameters {v} and {w} such that: 50—0 

1) Global (Finite element method) 

2) Local (Monte Carlo) 

We use the velocity Monte Carlo approach to find the set of variational parameters. 

(F. CtenJ>hysicaA282, 339, 2000) 

After each VMC step using {v} and {<»} update {r} and {(j)} using a simple forward integration 


Argom'ie Ncsioncd Sakaratary 
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DIMENSIONAL ANALYSIS; CHARACTERISTIC LENGTH AND TIME 
SCALES IN THE PRESENCE OF BOTH 
GB MIGRATION AND GRAIN ROTATION 


Grain growth due to either GB migration or grain rotation alone exhibits no physical 
length scale. However the presence of both GB migration and grain rotation introduces a 
physical length scale Rc This is determined solely by the material parameters and 
temperature 


Dimensional analysis; characteristic length and time scales 
in the presence of both GB migration and grain rotas ion 

Grain growth due to either GB migration or grain rotation alone exhibits no physical length 
scale. However the presence of both GB migration and grain rotation introduces a physical 
length scale R c This is determined solely by the material parameters and temperature 
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physical length scale 
(R c ~ 3.2 mn for Pd) 





Figure 14 


315 




TOPOLOGICAL CHANGES 


To account for the discontinuities occurring during the evolution of the 
microstructure certain topological changes has to be implemented. In addition to the well 
known neighbor switching (Tl) and grain disappearance (T2) events we have introduced 
an additional topological change: the grain disappearance by grain coalescence. 



Figure 15 
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GRAIN GROWTH IN THE PRESENCE OF BOTH GB MIGRATION AND 
GRAIN ROTATION. COMPARISON BETWEEN ATOMISTIC AND 
MESOSCALE SIMULATION OF A 25 GRAIN SYSTEM 


Here are shown side by side a few snapshots of the evolving microstructure. These 
are obtained from MD and mesoscale simulations. The parameters listed bellow and 
characterizing the mesoscale model were extracted from MD simulations: 

Ymax= 1.01 J m' 1 ; m max = 9.27xl0' 8 m 4 /(Js) ; Ro=14nm;R c = 3.2 nm; r| =(Rc/Ro ) 2 
= 0.05 

initial microstructure, t=0.0 ns 


Grain growth in the presence of both GE migration and grain rotation. 
Comparison between atomistic and mesoscale simulation of a 25 grain system. 

Here are shown side by side a few snapshots of the evolving micro structure. These are 
obtained from MD and mesoscale simulations. The parameters listed bellow and 
characterizing the mesoscale model were extracted from MD simulations : 

y max = 1.01 Jm- 1 ; m TTiax = 9.27x1 0 8 m 4 /(Js) ; R„= 14 nm; R c =3.2nm; n =(R c /R o) 2 = 0 05 
mitial mia ©structure, t-0.0 i-s 
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ATOMISTIC AND MESCOPIC 



Figure 17 
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VALIDATION OF GRAIN GROWTH SIMULATIONS BASED ON 
VARIATIONAL FUNCTIONAL FOR THE DISSIPATED POWER. 


Grain growth study in the presence of GB migration only: 

1. Isotropic GBs 

2. Anisotropic GBs 


Validation of grain growth simulations based on 
variational functional for the dissipated power. 


Grain growth study in the presence of GB migration only: 
1. Isotropic GBs 
2. Anisotropic GBs 





Figure 18 


319 





GROWTH LAW IN THE PRESENCE OF GB MIGRATION ONLY 


(D. Moldovan et al., Phil. Mag. A, 2002 (in press)) 

The kinetics of the growth process is characterized by the growth law. The average 
grain area is plotted as a function of time for both the isotropic and anisotropic systems. 
After some initial time the two systems reach the scaling regime, showing linearity in the 
long-time behaviour which is characterized by a grain-growth exponent n, that is, A(t) ~ 
t n . Specifically for the isotropic system the mean grain area increases linearly with time 
(i.e. n=1.0) while for the anisotropic system we obtain n=0.71. 


Growth law in the presence of GE migration only 

(D. Moldovan et al., Phil. Mag. A, 2002 (in press)) 

The kinetics of the growth process is characterized by the growth law. The 
average grain area is plotted as a function of time for both the isotropic and 
anisotropic systems. After some initial time the two systems reach the 
scaling regime, showing linearity in the long-time behaviour which is 
characterized by a grain-growth exponent n, that is, A(t) ~ t n . Specifically for 
the isotropic system the mean grain area increases linearly with time (i.e. 
n=1.0) while for the anisotropic system we obtain n=0.71. 
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GRAIN GROWTH CHANGES THE FRACTION 
OF LOW-ANGLE AND HIGH-ANGLE GBS 


The fraction of low-angle and high-angle GBs is being modified during growth. This 
is best characterized by the misorientation-angle distribution function, which basically 
gives the fraction of the total number of GBs with a given misorientation. The increase of 
the fraction of small-angle (low-energy, low-mobility) boundaries in the system explains 
the lower growth-law exponent for the anisotropic system compared with the isotropic 
case. This is mainly because the average driving force, proportional to <y(t)>/<R(t)>, 
decreases gradually with time. 


Grain growth changes the fraction of low-angle and high- 
angle GBs 

The fraction of low-angle and high-angle GBs is being modified 
during growth. This is best characterized by the misorientation-angle 
distribution function, which basically gives the fraction of the total 
number of GBs with a given misorientation. The increase of the 
fraction of small-angle (low-eneigy, low-mobility) boundaries in the 
system explains the lower growth-law exponent for the anisotropic 
system compared with the isotropic case. This is mainly because the 
average driving force, proportional to <y(t)>/<R(t)>, decreases 
gradually with time. 
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GRAIN-SIZE DISTRIBUTION FUNCTIONS 


The grain-size distribution functions for the isotropic and anisotropic systems 
averaged over a few configurations in the scaling regime are compared with the fitted 
log-normal distribution and the theoretical distribution functions of Hillert and Louat. It 
can be seen that both distributions obtained from simulations the cut-off value of the 
reduced grain size is about x 0 = 2.2. Moreover the shapes and the positions of the peaks 
are nearly the same for both the isotropic and the anisotropic systems. At small grain 
sizes they follow closely the Hillert’s distribution whereas at large grain sizes they 
decrease less steeply and come close to Louat’s distribution. 


Grain-size distribution functions 

The grain-size distribution functions for the isotropic and anisotropic 
systems averaged over a few configurations in the scaling regime are 
compared with the fitted log-normal distribution and the theoretical 
distribution functions of Hillert and Louat. It can be seen that both 
distributions obtained from simulations the cut-off value of the 
reduced grain size is about x c = 2.2. Moreover the shapes and the 
positions of the peaks are nearly the same for both the isotropic and 
the anisotropic systems. At small grain sizes they follow closely the 
Hillert’s distribution whereas at large grain sizes they decrease less 
steeply and come close to Louat’s distribution. 
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THE EFFECT OF GB DISCRETISATION ON THE TRIPLE-POINT ANGLE 

DISTRIBUTION FUNCTION 


For an isotropic system it is well known that the triple-point angles are all equal to 
120°, and this value is predicted by the Herring relation(given at the bottom of this page). 
However in our simulations with discretized GBs the adherence to Herring relation 
depends on the degree of discretization of each GB. The ability of each GB to assume a 
curved configuration depends on the number of boundary nodes along each boundary. 
The effect of the number of boundary nodes on the triple-point equilibrium angles can be 
understood best by plotting the triple-point distribution for systems with various numbers 
of boundary nodes. For the isotropic case the deviation from 120° is a measure of the 
effect of the degree of boundary discretization on the triple point-equilibrium condition. 
According to the Herring relation in the anisotropic system the triple-point angles may be 
substantially different from 120°. However, the extent to which the Herring relation is 
obeyed can be established by plotting the triple-point angle-deviation distribution 
function. 


The effect of GB discretisation on the triple-point angle distribution 
function 


For an isotropic system it is well known that the triple-point angles are all equal to 
120°, and this value is predicted by the Herring relation(given at the bottom of this 
page). However in our simulations with discretized GBs the adherence to Herring 
relation depends on the degree of discretization of each GB. The ability of each GB 
to assume a curved configuration depends on the number of boundary nodes along 
each boundary. The effect of the number of boundary nodes on the triple-point 
equilibrium angles can be understood best by plotting the triple-point distribution for 
systems with various numbers of boundary nodes. For the isotropic case the 
deviation from 120° is a measure of the effect of the degree of boundary 
discretization on the triple point-equilibrium condition. According to the Herring 
relation in the anisotropic system the triple-point angles may be substantially 
different from 120°. However, the extent to which the Herring relation is obeyed can 
be established by plotting the triple-point angle-deviation distribution function. 
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GROWTH LAW IN THE PRESENCE OF GRAIN ROTATION ONLY 


The effect of the grain rotation on the kinetics of grain growth are best revealed by 
investigating the time evolution of the average grain area. To emphasize the differences 
in growth the results obtained for r| = oo (grain rotation only) are analyzed. Our 
simulations show that the growth exponent n in the presence of grain rotation only is 
smaller than in the presence of GB migration only. Moreover the growth exponent 
depends on the actual grain size dependence of the grain rotational mobility considered in 
the simulations 


Growth law in the presence of grain rotation only 

The effect of the grain rotation on the kinetics of grain growth are best 
revealed by investigating the time evolution of the average grain area. To 
emphasize the differences in growth the results obtained for r|= oo (grain 
rotation only) are analyzed. Our simulations show that the growth exponent 
n in the presence of grain rotation only is smaller than in the presence of 
GB migration only. Moreover the growth exponent depends on the actual 
grain size dependence of the grain rotational mobility considered in the 
simulations 



M ( R ) — ^ p Grain rotational mobility 


Two accommodation mechanisms for grain rotation are considered: 
p=5 (GB diffusion) and p=4 (volume diffusion) 





Figure 23 


324 




GRAIN-COALESCENCE CHANGES THE FRACTION 
OF LOW-ANGLE AND HIGH-ANGLE GBS 


The fraction of low-angle and high-angle GBs is being modified during growth by 
grain rotation grain coalescence only. In particular one can see that there is significant 
decrease of the low-angle GBs fraction. 


Grain-coalescence changes the fraction of low-angle and high- 
angle GBs 

The fraction of low-angle and high-angle GBs is being modified during 
growth by grain rotation grain coalescence only. In particular one can see 
that there is significant decrease of the low-angle GBs fraction. 
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GRAIN GROWTH IN THE PRESENCE OF 
BOTH GB MIGRATION AND GRAIN ROTATION. 



Figure 25 
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GROWTH LAW IN THE PRESENCE OF 
BOTH GB MIGRATION AND GRAIN ROTATION (P = 5 MODEL) 


The effect of the grain rotation on the kinetics of grain growth are best revealed by 
investigating the time evolution of the average grain area. As shown in the figures bellow 
the value of the growth exponent, n, depends on the parameter h. Specifically we find that 
for r| <10, the growth exponent is practically equal to n=0.71, which is the same as the 
value for the system in the presence of GB-migration only. Similarly, for r|>3000 we 
find n=0.5, i.e., the same as the value obtained in the presence of growth by grain 
rotation-coalescence only. In the range 10 < rj <3000 , n decrease smoothly from n=0.71 
to n=0.5. This is the crossover regime where both growth mechanisms have comparable 
contribution to the growth process. 


Growth law In the presence of both GB migration and grain rotation 
(p = 5 model) 

The effect of the grain rotation on the kinetics of grain growth are best revealed 
by investigating the time evolution of the average grain area. As shown in the 
figures bellow the value of the growth exponent, n, depends on the parameter r\. 
Specifically we find that for r\ <10, the growth exponent is practically equal to 
n=0.71, which is the same as the value for the system in the presence of GB- 
migration only. Similarly, for ri > 3000 we find n=0.5, i.e., the same as the value 
obtained in the presence of growth by grain rotation-coalescence only. In the 
range 10 < r| <3000 , n decrease smoothly from n=0.71 to n=0.5. This is the 
crossover regime where both growth mechanisms have comparable contribution 
to the growth process. 
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MISORIENTATION DISTRIBUTION FUNCTION IN THE PRESENCE OF 
BOTH GB MIGRATION AND GRAIN ROTATION 


The fraction of low and high-angle GBs in the system evolves quite different when 
both GB migration and grain rotation are present. As previously shown during the growth 
by grain rotation-coalescence alone, the fraction of low-angle decreases significantly, 
while for curvature-driven GB migration alone the opposite trend is observed. In the 
crossover regime when migration and rotation have comparable contributions to the 
growth both trends are present. For example figure below shows for r|=250 that initially 
when the grains are small, rotation-coalescence dominates and the fraction of low-angle 
GBs decrease. However as the average grain size increases the GB migration starts to 
dominate and the fraction of low-angle GBs increases. 


Miso dentation distribution function in the presence of both GB 
migration and grain rotation 

The fraction of low and high-angle GBs in the system evolves quite 
different when both GB migration and grain rotation are present. As 
previously shown during the growth by grain rotation-coalescence alone, 
the fraction of low-angle decreases significantly, while for curvature-driven 
GB migration alone the opposite trend is observed. In the crossover regime 
when migration and rotation have comparable contributions to the growth 
both trends are present. For example figure below shows for r) = 250 that 
initially when the grains are small, rotation-coalescence dominates and the 
fraction of low-angle GBs decrease. However as the average grain size 
increases the GB migration starts to dominate and the fraction of low-angle 
GBs increases. 
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GRAIN SIZE DISTRIBUTION FUNCTION 


The effect of the presence of both GB migration and grain rotation on the grain size 
distribution function is illustrated bellow by plotting them for systems with h=co (grain 
rotation only), r|=250 (migration+rotation) and r|=0 (GB migration only). For 
comparison, we also show the log-normal distribution function fitted to the distribution in 
the presence of GB migration only (r|=0). The shapes and peak positions are different for 
the three distribution functions. The distribution function for r|=oo is narrower than for the 
others two systems and has a higher peak value; notice also that the distribution function 
drops more quickly to zero at small grain sizes, showing an apparent cut-off at x o =0.2. 
The presence of GB migration leads to a widening of the grain-size distribution and a 
shift of the peak positions towards larger reduced grain sizes. 


Grain size distribution function 

The effect of the presence of both GB migration and grain rotation on the 
grain size distribution function is illustrated bellow by plotting them for 
systems with r|=co (grain rotation only), r|=250 (migration+rotation) and 
r|=0 (GB migration only). For comparison, we also show the log-normal 
distribution function fitted to the distribution in the presence of GB 
migration only (r|=0). The shapes and peak positions are different for the 
three distribution functions. The distribution function for r)=co is narrower 
than for the others two systems and has a higher peak value; notice also that 
the distribution function drops more quickly to zero at small grain sizes, 
showing an apparent cut-off at >^=0.2. The presence of GB migration leads 
to a widening of the grain-size distribution and a shift of the peak positions 
towards larger reduced grain sizes. 
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CONCLUSIONS 


• Our MD simulations show that at least at the nanoscale grain sizes there are two 
equally important grain-growth processes: GB migration and grain rotation. 

• We have developed a mesoscale model for grain growth which incorporates GB 
anisotropy and grain rotation. 

• The presence of both GB migration and grain rotation introduces a physical length 
scale Rc into the system. If the average grain size is smaller than Rc, grain-growth is 
grain-rotation dominated; by contrast if the average grain size is larger than Ro, grain 
growth is dominated by GB migration. 


CCMCyiSIONS; 


* Our KID simulations show that at least at the nanosc-ale grain 
sizes there are two equally important grain-growth processes: 
GB migration and grain rotation. 

* We have developed a mesoscale model for grain growth which 
incotporates GB anisotropy and grain rotation. 

* The presence of both GB migration and grain rotation 
introduces a physical length scale R,, into the system. If the 
average grain size is smaller than R c , grain-growth is grain- 
rotation dominated; by contrast if the average grain size is larger 
than R c . grain growth is dominated by GB migration. 
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Anisotropic Materials Properties: 
Functional Forms for Multiscale Modeling 


Nicholas Bailey 
Cornell University 
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INTRODUCTION 


This talk is about a particular aspect of the "parameter passing" method of multiscale 
modeling, namely the mathematical representation of the parameters in question. 
Typically the information to be passed is more than a few numbers; it’s functions of 
several variables: orientation, stress etc. From our experience in three different contexts 
we have abstracted an approach to designing functional forms which is rooted in basic 
physics principles. The examples are (1) the stress dependence of Peierls barriers of 
dislocations (2) the orientation dependence of atomistically computed surface energies 
and (3) the orientation dependence of the etching rate of silicon by KOH. 


Anisotropic Materials Propsrtiss: Functional norms 



for Multiscale Modeling 




Nicholas Bailey Thierry Cretegny 

Markus Rauscher James Sethna 
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MULTISCALE MODELING 


In multiscale modeling we have atomistic simulations dealing with the smallest scales 
and some combination of continuum models and defect dynamics at larger scales. Going 
from atomistic-based modeling to defect based modeling means changing the 
mathematical description of the material; defects are primary entities rather than 
emergent entities. We reduce the number of degrees of freedom, but also in a sense 
increase the precision of the description, (e.g. by requiring a defect to have an exact 
location). This is an artifact of the change of description but it means that care must be 
taken in the formulation of defect rules. The rules are those of defect dynamics: e.g. 
mobilities as a function of all continuum parameters. The parameters come first from 
geometry: There are fewer entities, but it takes more numbers to specify their geometry 
(five numbers for a grain boundary). Then there is the interaction of a defect with its 
environment: effect fields at the position of the defect (stress, electric field, concentration 
of hydrogen etc) but also the contribution to these fields elsewhere in the material from 
the defect. 
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FUNCTIONAL FORMS PHILOSOPHY 


We have a philosophy for producing functional forms. The purpose is to create 
functional forms requiring the fewest fitting parameters, which can be fit with the fewest 
data points, given that there are typically many independent variables. There are five 
ingredients of the philosophy. Symmetry that is known to be present in the data, such as 
crystal symmetry (e.g. cubic), must be built into the functional forms. The same goes 
with singularities which are known to be in the data, such as cusps— a naive approach to 
fitting data would be to use an analytic expansion, but that would take many terms to 
reproduce a singularity. It makes much more sense to include this right in the functional 
form (compare quarter point elements in fracture modeling). To actually suggest a 
functional form one looks for a very simple physical model of the defect. A model too 
simple to be of quantitative use can often suggest appropriate functional forms. The final 
ingredients are subtleties in definitions (features of the continuum description that are 
missing from the atomistic description) and extra physics associated with singularities. 
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DISLOCATION GLIDE 


My first example is the calculation of the Peierls barrier as a function of stress for 2D 
edge dislocations. This is connected to dislocation glide (in a 2D crystal, glide at low 
temperatures is given simply by thermal activation over the Peierls barrier). We use 
flexible boundary conditions (boundary atom positions are given by linear elasticity 
including multipole terms) and the Nudged Elastic Band method to compute the barrier 
for a range of applied stresses. The symmetry is the periodicity of the lattice, which 
implies a relation between the forward and backward barriers (positive and negative shear 
stress). The singularity is the vanishing of the barrier at the Peierls stress as a 3/2 power 
law (saddle-node bifurcation). See next slide for the functional form. The subtlety is the 
definition of the center of the dislocation, which is connected with the multipole 
coefficients (a change of the center can be absorbed into a change in multipole 
coefficients, so there is a choice, like choosing a gauge in electromagnetism). The extra 
physics is that as the barrier vanishes at large stress, double and multiple jumps become 
important (before continuous sliding happens). 



□ Symmetry: Eg( a^) Eg(a- xy ) — I 

□ Singularity: Saddle-Node Bifurcation => X — — e -j- X 2 ■-■■■■> Eg •--- e 3 / 2 , , . 

□ Simple modal: 1 D sinusoidal potential E(x, a xy } — A sin(2rra.'/fe) ka^z 

□ Subtleties: Multipoles/position of dislocation ambiguous 

□ Extra Physics: Double (etc) jumps near depinning ( E % = 0) 


Figure 4 
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FUNCTIONAL FORM FOR BARRIER DATA 


Here is the functional form we are using. We consider a one-dimensional potential, 
with sinusoidal and linear terms, for which we can calculate the barrier exactly without 
the analytic factor h it is exact for the sinusoidal potential. Putting in the h allows fitting 
to real data. Choosing h in this way maintains the symmetry and the correct form of the 
power law singularity. Here is one possible choice of h; others can be imagined. Ideally 
only one or two terms in the h series are needed. The figure on the left shows a two 
parameter fit to the dependence on shear stress; on the right is a contour plot of sigma-xx 
and sigma-yy. The dependence on the other components of stress is nearly linear so we 
include it by making the parameters of the functional form shown here depend on those 
components and expand to quadratic order (so the total number of parameters is three 
times the number in the single-variable function). 
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SURFACE PROPERTIES: NEAR EQUILIBRIUM 


My second example, due to Thierry Cretegny, is about the orientation dependence of 
surface energies (computed from an empirical potential though they could be from an ab 
initio calculation). This function will be a single function of two angular variables. The 
symmetry in this case is cubic, so rather than the whole unit sphere we consider a 
triangular patch (figure on the left). Plotting energies along the indicated path gives the 
figure on the right. The dots are individually calculated energies; the line is a two- 
parameter fit to all the data in the triangular plot. Near high symmetry directions there are 
cusps in the surface energy (nearby surfaces are like steps on the high symmetry surface, 
the energy is proportional to the step density and hence the absolute value of the angle). 
To get a functional form we use a “broken bond model” which estimates the energy of a 
surface simply counting dangling bonds (ignoring relaxation of the atoms, which takes 
place in the actual computation). For an FCC lattice one needs only consider the cusp 
functions from the first and second neighbor bonds (110 and 100 directions). 


Surface Properties: Near equilibrium 

Thierry Cretegny: Surface Energies of Cu from simulation 



□ Symmetry: Cubic 


□ Singularities: Cusps at low-index surfaces 
ct(? 1) = A/ 11 o(n)+B/ 100 (n) 



/loo(ft) = 8(1*1 + \y\ + M) 


Broken-bond model-e cusp functions 



Figure 6 
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SURFACE ENERGIES: SUBTLETIES 


The subtlety in the case of the surface energies is that not all orientations yield a 
stable surface. Instead the surface tends to facet into two different surfaces. However for 
use in multiscale modeling, or say for making a Wulff plot to find the equilibrium crystal 
shape, it is desirable to have a surface energy for all orientations. To do this Thierry did a 
constrained minimization by choosing boundary conditions such that only one step fits. 
Extra physics is needed for example if studying growth of crystal surfaces in a continuum 
model. At the singularities in this data, the mobility of the surface vanishes, which 
corresponds to the fact that growth of the surface requires nonlinear effects coming from 
island nucleation etc. 


| J | ^ ^ Energies: Subtleties 

□ Subtlety: some surfaces are unstable! 

- for continuum model need all orientations 

- reconstructions, step bunching, faceting 

- prevent by only allowing one step within periodic boundary conditions 

□ Extra Physics: 

- Nonlinear mobility at cusps (facets have zero mobility in linear response) 

- Need nucleation of islands and pits, etc. 


Figure 7 
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SURFACE PROPERTIES: FAR FROM EQUILIBRIUM 


My last example is from Markus Rauscher and concerns modeling of experiments by 
Wind and Hines. They measure the etching rate of silicon by KOH within a particular 
crystallographic zone by fabricating a "pie" which is then cut into 1 degree wedges as 
shown. This creates a direct map of the etching rate within this zone (110 or 100). The 
modeling is to understand the faceting that occurs for certain orientations. His input was 
the measured rate, apart from the orientations where faceting occurs. The symmetry and 
singularities here are as for surface energies: silicon is cubic, and in a simple model the 
etching rate is proportional to the number of dangling bonds, i.e. a broken bond model 
again (now for the diamond lattice). As well as one cusp function (111) from the broken 
bond model he also added cubic harmonics (analogous to spherical harmonics but with 
the cubic symmetry built in). The middle figure shows the data from two different zones, 
along with the five parameter fit in those zones and the 111 zone. The figure on the right 
is a 3D plot of the fit. 


Anisotropic Etching Rates: Rauscher, Hinls, Wind " 



□ Symmetry, singularities, functional forms: As in equilibrium (broken bond 
models) 

n 5 parameters: 1 cusp function (< I'll >), 4 cubic harmonics (orders 0, 4, 6. 
8 ) 


Figure 8 
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ETCHING SUBTLETIES 


The subtlety here is very important; otherwise it might not make sense to use 
experimental data as input to a model whose outcome is to be compared to experiment. 
Because faceting appears in certain orientations, namely those near cusps as shown in the 
figure, the nominal measured rate there is an average over the facets and cannot be used. 
The data used is that corresponding to smooth surfaces and the cusps themselves. The 
model’s purpose is to understand the faceted regions and the transition from smooth to 
faceted, so there is indeed a separation between the input and the output of the model. 
The new physics needed here arises from the fact the PDE in the continuum model is not 
very well behaved and tends to form its own singularities; these are related to the faceting 
transition; when solved numerically without special treatment the solution blows up. This 
can be avoided by adding some diffusion (viscosity solution), but this leads to eventually 
flattening rather than facets. The problem was solved by including surface tension which 
comes into play near regions of high curvature. 


# ^ 1 1 1 *! 

x, v : : : : : : : : : : : x.x : x....x- x..x : x ; ..x- x. : : : x ..... x. x x. . ... ...:.: : x 



□ Subtlety 

- Fit only smooth surface data 

- Faceted surfaces average 

□ Extra physics: 

- Etching rate not enough at 
edges 

- Regularization-dependent: 

- Viscosity solution — s- 

rounding, eventual flattening 

- Surface tension solution 
(experimental anisotropic 
surface tension) -e- faceting 
as observed in experiment 


Figure 9 
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SUMMARY: GUIDING PRINCIPUES FOR USEFUU FUNCTIONAU FORMS 


To sum up, I want to reiterate that what we are presenting is more than a recipe, it is a 
philosophy, a way of thinking: symmetry, singularities, simple models, subtleties and 
extra physics. 


Summary: Guiding Principles 
for Useful Functional Forms 



& v fti rfii^trv 

to? f iililiv T 


Efficient 

Functional 



Forms 


subtleties 
new physics 

A V 



Figure 10 
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Properties of Inorganic Microporous Films 
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INTRODUCTION 


The focus of this presentation is on multiscale modeling in order to link processing, 
microscructure, and properties of materials 


Multi scale Modeling for Linking Growth, Microstructure, 
and Properties of Inorganic Microporous Films 

Dion G. Vlachos 

Department of Chemical Engineering and 
Center for Catalytic Science and Technology 

University of Delaware 
Newark, DE 19716 

vlachos@che . udel . edu 
www.che.udel.edu/vlachos 

Figure 1 
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Figure 2 


NUCLEATION AND GROWTH OF MATERIALS 


Overview of problems we study (clockwise): Growth mechanisms in chemical and 
physical vapor epitaxy, thin films of zeolites for separation and sensing, thin Pd films for 
hydrogen separation, and pattern formation by self-regulation routes 
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CHEMICAL REACTORS EXHIBIT MULTIPLE SCALES 


Schematic of multiple length scales in a deposition process. Three broad scales can be 
distinguished: The reactor scale, where the continuum approximation is valid, the 
mesoscale where the microscopic mechanisms control growth mode and crystal 
morphology, and the quantum scale. Quantum scale controls intermolecular forces which 
in turn affect nucleation and the kinetics of diffusion and surface reactions. Phenomena at 
each scale are coupled in a two-way manner. From a controls point of view, solving at 
least the reactor and mesoscale models simultaneously is necessary. 



Figure 3 
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OVERALL MULTISCALE MODELING APPROACH 


Overall approach we undertake. First, three levels of models are bridged together: 
chemistry from semi-empirical methods (especially for large reaction mechanisms) or 
DFT, kinetic MC for the surface, and continuum models for the reactor. Through semi- 
empirical parametrization or lookup tables, we decouple the quantum and mesoscale 
simulations. Meso- and macro-scale models are fully coupled. Methods for achieving this 
are described in Raimondeau and Vlachos, Chem. Eng. J. (accepted). 

In a number of problems, kinetic MC is inadequate to capture realistic length scales. 
Mathematically rigorous coarse-grained models are then derived from the master 
equation in order to replace MC in multiscale frameworks. We present an example of this 
approach in this talk. 

Finally, control of microstructure of materials demands reduction of these complex 
models to a small number of ODEs. An example is discussed in this talk. 




• Efficient direct numerical simulations 

of coupled reactor scale-microscopic scale 
models are developed 


theory and 
macroscopic laws of pattern 
evolution are derived 


* Reduced models are developed for efficient 
design and on-line control 


Reactor scale model: 
Quantum mechanics Fluids, heat, species, 
detailed kinetics 

\ // 

Monte Carlo 


Master equation 


'v&ubsc stoic scale: 

Pattern evoi ution laws 
via transport and 
thermo properties 


Mesoscopic scale: 

Ix>cal mean field theory 


• These models have great potential for rational 
design of novel materials (e.g., defect free zeolite 
membranes, heterostructures with large 
mismatch) 


Figure 4 
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FULL COUPLING BETWEEN FUILD (MACRO) AND SURFACE (MESO) 

SCALES 


The first example discussed is a simple model system that illustrates the coupling of 
meso and macro scales in a vertical reactor. Spatially average rates are computed at each 
time step of MC and fed into the boundary conditions of the continuum model, which in 
turn provides the concentrations at the gas-solid interface needed in the MC. 



Gas-phase continuum model / Finite differences 



Solid-on-solid lattice Surface parameters 

First-nearest neighbors interactions Rate constants, surface T,... 

Adsorption-Desorption-Migration 

Vlachos, Appl. Phys. Lett. 74, 2797 (1999) 

Figure 5 
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MACROSCOPIC PARAMETERS CONTROL MICROSTRUCTURE 


Example of coupled, multiscale growth simulations. These graphs show the effect of 
substrate temperature on growth rate and microroughness at three strain rates 
(proportional to flow rates). At high temperatures step flow is observed, whereas at lower 
temperatures, 2D nucleation dominates giving rise to increased roughness. These 
simulations provide for the first time insight into the role of macroscopic, experimentally 
accessible quantities, such as temperature and flow rate, into the transition from one 
growth mode to the other (demarcated by open circles). 
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PROPER ORTHOGONAL DECOMPOSITION (POD) 


Model reduction is essential from an engineering perspective in order to develop 
efficient models that can be used for on-line control. The POD method is used for fluid 
phase and surface data reduction obtained from hybrid, multiscale simulations. 


Proper orthogonal decomposition (POD) 


• POD is used to reduce temporal or spatiotemporal data of complex 
dynamic models to low-dimensional systems (Sirovich et al., 1987) 

• {u (k) :k=l , . . . ,M} is a set of M snapshots from M1H algorithms 

• N is the number of modes necessary to capture the dynamics of the 
system based on the eigenvalues 


Time correlation matrix (< 

=fy D rf%x)tP(x)dx, U = 1 

Eigenvectors of (C) 


Eigenfunctions 

M 

<*>„(*) Z4 /n "' ;, (-0 n=\,...,M\ 

k 1 


SVD -> a n (t) 


Reduced system 

! N 

K t > x )=l^ a n( t )k(x) + U(x)\ 
n- 1 I 


Figure 7 


351 











RESULTS OF MODEL REDUCTION 


Example of model reduction. A step change in the precursor concentration is 
introduced, and the dynamics of the system is monitored. The reduced model (~6 ODEs) 
does an excellent job for the fluid phase. Surface morphology is also well captured, but 
some smearing does occur. DNS stands for direct numerical simulations and POD 
represents the reduced model. 






3 4 

711116 [s] 


" 1 15 

Reduction of fluid phase data 

- | 

in a stagnation flow geometry, 

- 5 a 
& 

C Q 3 

upon a step change in the 

--5 £• 

precursor’s mole fraction, is 

" • 10 s. 

excellent. 

■■15 I 

Raimondeau and Vlachos, 

J-2Q 

J. Comp, Phys. 160, 564 (1 


Reproducing film microroiighness with a low-dimensional model 
is important for on-line control. 




1 monolayer 


DNS 


T=1100 

a=l s' 1 
t=2.7 s 


Figure 8 
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RESULTS OF MODEL REDUCTION (CONTINUED) 


Example of model reduction at low temperatures. Now, 2D nucleation prevails. The 
reduction is again reasonable. The bottom right graph shows the average microroughness 
as a function of time. Upon introduction of a step change in the precursor, a new quasi- 
steady state is attained. The reduced model captures this information reasonably well, 
despite the fact that it was not directly trained for this task. Such data can be collected 
experimentally by STM or AFM. The performance of the reduced model indicates that it 
would be possible to use a STM or AFM tip as a sensor and reduced models derived from 
first principles, multiscale simulations, for on-line control of microstructure. 


Results of model reduction (Cont.) 


Film morphology from DNS and POD at t=2.7 s 




T=600 K 
a=100 s-‘ 

iocu 

98- 
96 
94 

1 monolayer ~ 0.5 nm 


DNS 


98 


Reduced models can reasonably 
capture film mierorotighness 

during dynamic changes in 
operating conditions. 
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<-v — 

20 
15 




Step change in piecnisoi mole fraction 

H V V-i 




DNS 
- - - POD 


Figure 9 
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MULTISCALE MODELING OF ZEOLITE MEMBRANES 


Introduction to zeolites and their applications as another example of multiscale 
simulations. Summary of three-step synthesis of zeolite films developed by M. Tsapatsis 
and co-workers. Controlling properties of membranes demands an understanding of the 
fundamentals of growth, which are discussed next. 



• Zeolites are inorganic, crystalline 
microporous materials used 
extensively in separations, reactions, 
optico-electronic materials growth, etc. 

• Multistep synthesis is followed 

• Hydrothermal synthesis of <0. 1 pm particles 

• Deposition of these particles on a substrate 

• Secondary growth of interparticle voids 

• Optimizing membrane microstructure 
requires understanding of: 

• Colloidal interactions and growth 
mechanisms 

• Particle morphology-processing conditions 
relations 


Formation of defects during the evolutionary precursor film 
roblem of internarticle growth 


gtrx 






SqV 




z z \ X 

•Yyr 
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SILICALITE-1 SYNTHESIS FROM CLEAR TPAOH-TEOS-H 2 0 SOLUTIONS 


Synthesis of zeolite nanoparticles takes place in aqueous solution, by adding a Si0 2 
source and a template (a structure directing agent). Hydrolysis and condensation reactions 
lead to a plethora of siliceous species and nanoparticles seen by a variety of experimental 
techniques indicated. The growth unit that contributes to growth has been, however, 
speculative till recently. 
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MODELING OF GROWTH 


To understand the growth and isolate it from nucleation, DLS experiments have been 
carried out for seeded growth and modeled with a transport model based on the chemical 
potential formulation. The chemical potential considers the Brownian motion of growth 
units and interactions between particles with the DLVO potential. Good agreement is 
found. Model parameters have been obtained from independent experiments in this case. 
A primary outcome of this work is that the growth unit appears to be nanoparticles of ~3 
nm in diameter. 


Modeling of growth 


DLS data are modeled at the nm scale 

250 r 
200 ‘ 




^Soor^oC oOOOOuOUo UuUuOOO 


2 4 6 

Time [hours] 

DLVO interactions play a key role in 
controlling growth rale and activation 
energies 

Nikolakis ctal., Chan. Mat. 12,845 
( 2000 ) 


10 


The model describes 
colloidal interactions during 
transport of nanoparticles to 
a zeol ite crystal usi ng tl: 

DL VO potential 


Growth rate is well 
The growth unit is ~ 3 nm 


Model parameters 
from 


Figure 12 
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AFM MEASUREMENTS OF FORCES VALIDATE DLVO POTENTIAL 


AFM experiments between a glass sphere and a zeolite film have been used to 
validate the DLVO potential used in the growth model. The force can be described 
reasonably by the DLVO potential. 


AFM measurements offerees validate 
DLV O potenti al 



Figure 13 


357 





TEM INDICATES THE EXISTENCE OF SUBCOLLOIDAL PARTICLES 


High resolution TEM indicates the existence of nanoparticles predicted by the model. 
Models and experiments have been used in a symbiotic manner for model validation. 



Figure 14 
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MODELING AT VARIOUS SCALES 


Armed with a growth model that describes growth at the nanometer scale, we are 
currently exploring ways to predict particle anisotropy and particle morphology evolution 
as a function of processing conditions. It turns that particle morphology strongly affects 
membrane microstructure. 


Modeling at ¥arious scales 


DLS data are modeled at the nm scale 



DLVO interactions play a key role in 
controlling growth rale and activation 
energies 

Nikolakis etal. Chan. Mat. 12,845 
( 2000 ) 


Anisotropy of crystals at 1. tun 


Membrane Shape of crystals 

orientation synthesized in solution 



Preferred orientation in zeolite 
membranes depends on the growth 
rates of crystallographic planes, 

j Bonilla etal., Micro. Meso. Mat. 42, 191 

Hi 


Figure 15 
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ZEOLITE MEMBRANE GROWTH IS MODELED 


Fundamental growth models can be used to study the evolution of thin films from an 
initial coating of particles. The left graph shows 2D simulations of 10,000 particles as a 
function of time or height. The right graph shows a cross section SEM of a silicalite-1 
zeolite film. The model predicts the columnar structure of films. 




5 10 15 20 25 30 

(Hm) 

2D Modeling 


Cross Section SEM 


• Columnar grain boundary structure 

• Average grain size increases with increasing film thickness 


Bonilla etal.. Micro. Meso. Mat 42, 191 (2001) 


Figure 16 
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DEFECTS CAN BE PREDICTED AND CONTROLLED 


Defects can also be predicted that happen near the zeolite/substrate interface. 
Simulations show that the most critical factor controlling defects is the growth factor, a, 
which is the ratio of velocities of different particle facets. Under more isotropic growth 
conditions, defects are minimized. Thus, it is important to link processing conditions with 
particle morphology. 


Defects can be predicted and controlled 

2.5 

col.35 

2.0 

1.5 

I 

X 1.0 

0.5 
0.0 

33 33.5 34 34.5 35 35.5 36 

x (pm) 

Minimize defects by close packing of seed coating and control of initial roughness 

Figure 17 
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XRD POLE FIGURE RESULTS 


Experiments for the primary orientation of these films indicate that either 2 or 1 
peaks, depending on growth conditions. 



Figure 18 
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MEMBRANE ORIENTATION IS AFFECTED BY PARTICLE 
MORPHOLOGY AND GROWTH CONDITIONS 


The model can explain the formation of single or double peaks. The key again is the 
growth factor. These simulations provide insight into how to control the orientation of 
polycrystalline films that can have a significant impact on performance. 



Figure 19 
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DIFFUSION/REACTION THROUGH MEMBRANES 


Confocal micrograph of zeolite films indicates the presence of grain boundaries. 
Prediction of transport of molecules through these films is a challenging, multiscale 
problem that direct molecular simulations cannot handle. Mesoscopic theory emerges as 
an invaluable tool in bridging the length and time scales problem in predicting properties 
of polycrystalline films. 



gapill 




membrane/support 
interface (A1 2 0 3 ) 


Laser confocal micrograph of MFI zeolite membranes, 

Bonilla el al., J. Membrane Sci. 182, 103 (2001) 


Polycrystalline medium 
Intermolecular forces between adsorbates 
Non-equilibrium transport 
Scales of several microns 


Atomistic apprc 
is not adequate 


pic theory 


Figure 20 
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MESOSCOPIC THEORY OF PATTERN EVOLUTION 


Mesoscopic theory is developed from the microscopic master equation by non 
equilibrium statistical mechanics. This mesoscopic theory captures the details of 
microscopic mechanisms and is exact when the interaction potential is infinite ranged. 
Thermodynamic and transport macroscopic equations (e.g., velocity of moving 
interfaces, mobility, critical nucleus size, surface tension) can also be derived by further 
coarse-graining using homogenization theory and asymptotic expansion of traveling 
waves. 


Mesoscopic theory of pattern evolution 


Microscopic models (master equation) 


c scale coverage 
ge numbers 


Law of large numbers 



Mesoscopic theory 


d t u = V • \p\Vu - J3u(l - u)V(J *w)]] 


Asymptotic expansion using traveling 
waves and homogenization techniques 


l ickian term 


“Interaction” tenn 


Macroscopic laws of evolution 

- Normal velocity v = a™ + ( l\ 

- Mobility (Kubo-Green formulae) v=P 

- Surface tension 


J -“M 1-9) f 


Arrhenius dynamics 

v(x —>y)=d exp (~/3U(x)) 

D=d e- Pu 


a= \L\ J { r ' Mt +r ‘ ■«)*(£)(«•'■' i dr ' d t Katsoulakis and Vlachos, Phys. 

Rev. Lett. 84, 1511 


Figure 21 
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VALIDATION OF MESCOSCOPIC THEORY 


Numerical simulations in thin, single crystal films show excellent agreement between 
gradient (no equilibrium) MC simulations and mesoscopic simulations for various 
strengths of the interaction potential for relatively short ranged potentials. This result is a 
consequence of a Large Deviation Principle that these equations satisfy. 


Validation of mesoscopic theory 


1 -D non equilibrium problem 



Vlachos and Katsoulakis, 

Phys. Rev. Lett. 85, 3898 (2000) 



• Mesoscopic theory is compared to MC simulations with excellent: results 
for different interatomic potentials. 

• Sharp gradients in concentration occur for strong attractive interactions 
or low temperatures within the spinodal regime. 

• Effective diffusivities are strongly affected by mtermolecular forces. 

• Mesoscopic solutions are obtained extremely efficiently. 


Figure 22 
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BREAK OF MEMBRANE SYMMETRY 


Extension of mesoscopic theory from single crystal to polycrystalline materials can 
be achieved using homogenization theory. Simulations depicted here indicate that both 
the concentration and location of defects are important factors in controlling permeation 
properties. For example, the performance of a membrane depends on which side is 
exposed to the high pressure (feed side) when intermolecular forces dominate transport of 
species through a membrane. 


Break of membrane symme try 

Defects pattern: 


• In the presence of intermolecular 
forces, defects are more important 
near the high pressure side of a 
membrane 

•The flux varies in a nonlinear 
way with defect density 

Less defects 


More defects 


Figure 23 
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• Mesoscopic equations have been validated by comparison of 
MC to mesoscopic theory (FD and Spectral solutions) for a 
variety of potentials, microscopic dynamics, steady state and 
transients 

• The excellent performance of mesoscopic theory is a result of a 
Large Deviation Principle 

• Mesoscopic equations can be used to describe the behavior of 
heterogeneous systems (polycrystalline membranes, catalysts, 
composite materials) in the presence of intermolecular forces by 
extending homogenization techniques 

• Mesoscopic theory is an attractive alternative to atomistic 
simulations for larger scales and heterogeneous media 


Figure 24 


Summary 


• We show the feasibility of applying multiscale hybrid 
(quantum/atomistic/continuum) models to homogeneous- 
heterogeneous reactors 

- Surface chemistry developed assures thermodynamic consistency and 
coverage dependence of kinetic parameters. 

- Surface reaction mechanisms with predictive capabilities are developed. 

• Reduction of multiscale models is feasible and enables design and 
control of materials’ microstructure at the atomic level 

• Mesoscopic theory is a viable substitute to traditional molecular 
simulations for larger time and length scales 

- The theory compares well to Monte Carlo simulations in ID and 2D even for 
short range interactions. 

- Solution of mesoscopic theories are obtained in seconds. 


Figure 25 
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TALK OUTLINE AND COAUTHORS 


In the first part of this presentation, a combination of atomic-scale structural 
modeling, density functional calculations and mesoscale theory will be presented that 
yields predictions of the mechanical properties of polycrystalline materials. This is an 
example of hierarchical multiscale modeling in which results at one scale are used as 
input for methods appropriate at the next larger scale. This project is collaboration 
between Donald Brenner of North Carolina State University, Airat Nazarov and Alexei 
Romanov, of Russian Academy of Science, and Lin Yang of Lawrence Livermore 
National Laboratory. The Office of Naval Research provided support for this work. 

At the end of the presentation, a multiscale model of polymer nanocomposites which 
is another example of serial approach will be briefly outlined. The hierarchy of 
computational methods will include atomistic, mesoscopic and analytic levels. This work 
is a collaboration effort between Leonid Zhigilei, Dana Elzey and Ioannis Chasiotis of 
University of Virginia, Deepak Srivastava of NASA Ames Research Center and 
computational group of International Technology Center. 



Talk Outline & Coauthors 


Multiscale Modeling of Polycrystalline Materials 

Donald Brenner (NCSU) 

Airat Nazarov and Alexei Romanov (RAS, Russia) 

Lin Yang (LLNL) 


Multiscale Modeling of Carbon Nanotubes/ 
Polymer Composites (Introductory remarks) 

Leonid Zhigilei, Dana Elzey, Ioannis Chasiotis (UVA) 

Deepak Srivastava (NASA Ames Research Center ) 


Figure 1 
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MULTISCALE MODELING OF POLYCRYSTALLINE DIAMOND: 

ULTIMATE GOAL 


Accurate first-principles predictions of macroscopic-scale materials properties strictly 
from theory remain a long-standing goal of the materials modeling community. 
Regarding polycrystalline materials, it would be beneficial to accurately predict 
properties of their grain boundaries (GBs), since it is well recognized that many of the 
properties of polycrystals depend to a large extent on the properties of their internal 
interfaces and their junctions. GBs contribute to strength, ductility, toughness, corrosion 
resistance as well as thermal and electronic properties. In our research we concentrated 
mostly on prediction of mechanical properties of poly crystalline diamond according to an 
interest of ONR, however the developed approach is quite general. 



Macroscale mechanical 
properties - 


stability, toughness 


ft im 


Atomic scale - 
bonding, structure 
energy, stress 


Multiscale Modeling of Polycrystalline Diamond 

Ultimate Goal 


Real World: CVD Diamond 


Shechtman et.al. 


Strictly from theory; 
no experimental input! 


Figure 2 
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MULTISCALE MODELING OF POLYCRYSTALLINE DIAMOND 


While density functional methods are capable of accurately predicting properties of 
ideal structures and selected localized defects for many systems, they remain too 
computationally intensive to model extended defects typical of macroscopic-scale 
engineering materials. In contrast, continuum-based methods are ideal for accurately 
describing long-range stresses but are inappropriate for modeling atomic-scale chemical 
dynamics from which processes such as crack propagation originate. Intermediate 
between these two length scales is atomistic modeling using analytic interatomic 
potentials. While powerful, especially given current computing capabilities that allow 
very large systems to be modeled, this approach can yield results for which the range of 
validity is difficult to access due to possible inaccuracies in the interatomic potential. 



Continuum Models/ Finite Element Methods 

• Well developed for microstructure evolution 

• but ... missing atomic-level 'chemistry' of fracture 


Multiscale Modeling of P< 


How do we attack this problem? 

Direct First Principles Methods (e.g. density functional theory): 
•Accurate total energies 

•But ...limited to ~1000's of atoms on supercomputers! 


Analytic Potentials (e.g. EAM or bond order): 

•50,000 atoms on workstations/ 10 6 -10 8+ on parallel com puters 
•qualitative stress/energies and dynamics 
•but ... not always quantitatively accurate 


Figure 3 
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MULTISCALE MODELING OF POLYCRYSTALLINE DIAMOND 


That is why we used the following combined approach. In this scheme, the 
mesoscopic model gives a way to predict the energies and stresses for GBs in the entire 
misorientation range for a given tilt axes based only on a few reference structures of low- 
sigma GBs. These structures are small systems accessible to the first principles (FPs) 
density functional calculations. This leads to a completely first-principles estimate of the 
energy of tilt grain boundaries at any arbitrary tilt angle. Analytic potential calculations 
are used in the present scheme to prove the validity of dischnation approach to the 
covalent materials; generate starting configurations for the FPs calculations and to 
evaluate the extension of stress fields and consequently the sizes of computational cell for 
FPs calculations. In addition, large-scale molecular dynamics simulations with analytical 
potentials provide qualitative insight on dynamics of fracture initiation and propagation 
in polycrystalline materials. 



Classical 

Potentials 


Density Functional 
Theory 


initial structures 


system sizes 


Mesoscale Theory 

(Disclination/Structural Units Model) 


qualitative 
chemical 
dynamic; 


validate 

DSUM 


a few input 
energies 


Figure 4 


376 







ATOMISTIC SIMULATION METHODS 


The analytic potential energy function used in the present simulations is a many-body 
bond order potential developed by Donald Brenner that mimics bonding in carbon 
systems for a wide range of atomic hybridizations. The function reproduces a relatively 
large database of molecular properties of hydrocarbons as well as solid-state properties of 
carbon, including the in-plane lattice constant, cohesive energy and elastic properties of 
graphene sheets as well as cohesive energy and elastic properties of diamond. Density 
functional with local density approximation calculations had been carried out by Lin 
Yang of LLNL. 


Atomistic Simulation Methods 


, Inalvtic Bond Order Potential for Hydro-carbons 

- Second moment approximation to local DOS 

- Near neighbor pair sum modulated by analytic bond order 

- Empirical corrections for radicals, conjugation, dihedral rotation 

- Fit to extensive solid state and molecular data base 


First Principles Calculations 

Density functional local density approximation 

- Plane wave basis/60 Ry cut-off 

- Norm-conserving pseudopotentials 

- Effective k-point sampling 

- Relaxation of nuclear degrees freedom 
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DISCLINATION-STRUCTURAL UNITS MODEL 


Disclination-structural units model (DSUM) [1] combines the disclination description 
of GBs and structural imits (SUs) concept. A GB, being a rotational defect, is described 
more simply by disclinations, the liner rotational defects. A wedge disclination forms 
when one extracts or inserts a sector of the material from cylinder or disk. According to 
the disclination model, a general GB can be considered to consist of alternating areas of 
two special GBs. The jimctions of these areas are partial disclinations, and the boundary 
is represented as a disclination dipole wall (DDW). However, it was unknown how the 
period of DDWs scales with the geometrical period given by a coincidence site lattice. 
Later atomic calculations by Sutton and Vitec, which lead to the structural units model, 
demonstrated that the decomposition of general GB into sections of two low energy 
boundaries occurs within the geometrical period. All intermediate boundaries in a range 
delimited by two neighbor special boundaries are composed of the mixture of their SUs 
in a certain proportion and arrangement. The disclination formalism was subsequently 
combined with the structural units model to calculate grain boundary energies for metals 
[ 1 ]. 

[1] V. Gertsman, A. A. Nazarov, A.E. Romanov, et.al. Philos. Mag. A 59, 1113 
(1989) 
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GB ENERGY CALCULATIONS USING DSUM 


According to DSUM, at the lines along a GB where SUs of different types meet, the 
local misorientation angle alternates between 91 and 92 and therefore, these lines are 
partial wedge disclinations and the minority SUs can be considered as disclination 
dipoles. Because the DDW period is known from the coincidence site lattice parameters 
for the given boundary, its elastic energy is explicitly known. The energy of a GB is then 
given as a weighted sum of the energies of the individual SUs plus the elastic energy of 
the DDW and the energy of the disclination cores. This leads to a set of equations 
involving only the energies of the key single-SUs grain boundaries, bulk elastic 
properties, and a single unknown disclination core parameter. This parameter can be 
uniquely determined from the energy of a single GB containing a combination of 
structural units. The whole set of input parameters can be obtained from the first 
principles simulations leading to the accurate prediction of GBs energetics for the entire 
misorientation range for a given tilt axis. It should be noticed, that dislocation-structural 
units multiscale model for prediction of GB energetics had been developed by Wang and 
Vitek (1986), but it requires larger number of input parameters then DSUM. 
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VALIDATION OF THE MODEL: GB ENERGY CALCULATION 
WITH ANALYTICAL POTENTIAL 


To demonstrate accuracy of the DSUM approach for covalent materials, energies 
obtained directly from the bond-order potential for <001> symmetrical tilt GBs in 
diamond at different misorientation angles (indicated by the circles in the figure) were 
compared to values predicted by the DSUM (solid curve) [1], As input parameters for the 
DSUM the energies of two S=1 GBs, two S=5 GBs and S=29 GB from analytic potential 
calculations had been used (solid circles). Very good agreement can be seen between the 
energies predicted with DSUM and those calculated with the analytic potential. This 
demonstrates the validity of the disclination approach to covalent crystals. 

[l]O.A.Shenderova, D.W.Brenner , A.I.Nazarov, A.E. Romanov, L.Yang. Multiscale 
modeling approach for calculating grain boundaries energies from first principles', 
Phys.Rev.B, 57, R3181 (1998). 
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ACCURATE PREDICTION OF GBS ENERGIES 
USING FIRST PRINCIPUES SIMUUATIONS 


Below is GB energy versus misorientation curve obtained with input energies of 
special GBs and one intermediate GB calculated from the density functional approach. 
Energies of special GBs calculated from first principles are about 1.3-1. 5 times less those 
calculated with bond order potential due to electronic system relaxation. Necessary 
elastic parameters for diamond also had been calculated using density functional 
approach. So taking first principle calculations only for a few special GBs we obtain 
realistic GBs energies in the entire range of misorientation for a given tilt axis. 



Figure 9 
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USING DSUM FOR STRESS CALCULATIONS: COMPARISON WITH THE 
STRESS DISTRIBUTION CALCULATED USING THE ANALYTIC 

POTENTIAL 


In addition to GB energy, the disclination approach is able to provide stress 
distribution near GBs. Stresses near defected structures may also be calculated at the 
atomic level from interatomic forces and the result compared to the disclination model. 
The stress contribution attributed to individual atoms consist of “moments of the forces” 
between a central atom and its neighbors divided by the local atomic volume. Below such 
analysis is carried out for the S=149(7 10 0) GB, which has a relatively complicated 
stress field and therefore is a good test model [1], As it seen from the figure, small areas 
of compressive stress are inserted to more extended regions of tensile stress and vice 
versa. 

[1] O .A.Shenderova, D.W.Brenner, 'Atomistic Simulation of Grain Boundaries, 
Triple Junctions and Related Disclinations', Solid State Phenom., Trans Tech Publ., v.87, 
p.318 (2002) 
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CALCULATION OF STRESSES USING DSUM 


Figures below illustrate stresses in the vicinity of S=149 GB calculated using DSUM 
[1] and those calculated with the analytic bond-order potential. As it can be seen, very 
good agreement exist between the atomistic and mesoscopic approaches in the 
description of both shape of the stress distribution as well as stress magnitudes at 
distances from the interface more then about two interatomic distances. 

[1] R. Valiev, V. Vladimirov, et.al. “Stress fields of equilibrium and nonequilibrium 
grain boundaries”, Fiziko-technical inst. By Ioffe, Leningrad, 1989 (in Russian). 
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FIRST PRINCIPLE CLEAVAGE ENERGIES AND TOUGHNESS 
OF GBS IN DIAMOND 


GB energy plays a central role in various GBs properties, such as GB mobility and 
fracture, impurity segregation, GB diffusion, to name but a few. The simulation data for 
first principle GBs energies combined with the energies of free surfaces (also calculated 
from first principles) allow to evaluate cleavage fracture energy of the internal interface 
[1], The resulted curve is shown in the Figure. For comparison, the bulk ideal crystal 
cleavage energies are presented too. It is apparent from the Figure that special GBs 
possess higher cleavage energies relative to GBs in the nearby misorientation range. 
Cleavage energies of most <00 1> tilt GBs are about 60-75% of those for the ideal bulk 
crystals with the same orientation. Based on the cleavage energies, in turn, toughness of 
different GBs can be also evaluated. Results of calculations demonstrate good agreement 
with experimental data for CVD diamond fdms. 

[1] O.A. Shenderova, D.W.Brenner, A.Omeltchenko, X.Su, L.H.Yang. ‘Atomistic 
Modeling of Grain Boundaries Fracture in Diamond’, Phys.Rev . B 61 , 3877(2000). 
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EXTENSION OF DSUM TO MORE COMPLICATED GB STRUCTURES. 


GBs with a multiplicity of structures have been observed in Si and Ge. Some of them 
have been found to consist of zig-zag arrangements of structural units and others to 
contain dislocations having a screw component along the tilt axis. The <00 1> 
symmetrical tilt GBs commonly observed in materials with a diamond cubic lattice can 
be constructed in the entire range of misorientaton from three characteristic structural 
units: type A, the core of a pure edge dislocation, type B, the core of a mixed 45° 
dislocation, and units of perfect crystal. There exist two distinct groups of <00 1> tilt GBs 
that are delimited by the £=5/0=53.1° boundary. Grain boundaries with minimum 
energies in the misorientation range O0<<53. 1° consist of mixture of the structural units F 
and A s while GBs in the 53.1°0<<9O° misorientation range can consist of the F', A z and 
B s units. An intermediate GB can be composed of non-dissociated dislocations consisting 
of pairs of A z or B s units, or of separated A z and B s units (see Figure below). All types of 
structures are stable, although the energies are different. DSUM had been extended to 
treat these more complicated GBs structures. 
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EXTENSION OF DSUM TO MORE COMPLICATED GB STRUCTURES. 


In order to calculate GBs energies for groups of <00 1> GB models in the 
53.1°9<<90° misorientation range these structures were described in terms of flat and 
faceted disclination dipole walls and screw dislocation dipole walls (figures below), with 
the energies of these defects calculated from anisotropic elasticity theory. Details of the 
analysis are provided in [1], 

[1] A Nazarov, O.A. Shenderova, D.W.Brenner, ‘Elastic Models of <001> and <01 1> 
Tilt Grain Boundaries in Polycrystalline Diamond’, Phys.Rev.B 61 , 928(2000). 
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EXTENSION OF DSUM TO MORE COMPLICATED GB STRUCTURES. 


Resulted GBs for groups of metastable GBs predicted with DSUM are illustrated in 
the Figure below. Grain boundary energies for the <001> and <01 1> tilt axes for 
diamond given by atomistic modeling quantitatively agree with those from the DSUM for 
the most stable structures. For metastable structures, the DSUM predicts the correct 
energy ranking for different models at a given tilt angle, although the agreement between 
the results of the DSUM and atomistic model is not quantitative. The calculations 
demonstrate that an appropriate treatment of the elastic energy of complicated structures 
(e.g. facetted grain boundaries and grain boundaries with dislocations containing screw 
components) is needed to correctly predict energy rankings [1], 

[1] A Nazarov, O.A. Shenderova, D.W.Brenner, ‘Elastic Models of <001> and <01 1> 
Tilt Grain Boundaries in Polycrystalline Diamond’, Phys.Rev.B 61, 928(2000). 
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DSUM CALCULATIONS OF GBS ENERGIES FOR OTHER MATERIALS 


Examples below illustrate results of application of DSUM to <001> and <1 1 1> GBs 
in cooper [1], Analytic potential calculations of energies of delimiting GBs had been used 
as input parameters. Good agreement can be seen for results of simulations for a variety 
of GBs calculated with an analytic potential and those predicted by DSUM. Even if to use 
only S=1 delimiting boundaries (no atomistic simulations are required, interface energy is 
zero), still obtained curved posses very high predictive capability. 

[1] A Nazarov, O.A. Shenderova, D.W.Brenner, ‘On the disclinati on-structural unit 
model of grain boundaries’, Mat.Sci.Eng.A, 281,148(2000). 
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SUMMARY ON THE DISCUINATION-STRUCTURAU UNITS MODEU. 



Multiscale Modeling of Polycrystalline Materials 


Summary 


• Disclination -Structural Units Model (DSUM) allows 
accurate and internally validated first principles GB 
energies for the entire misorientation range 

• Energies of multiple GB structures can be calculated 
using DSUM 

• DSUM is valid for different types of materials 
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GRAIN BOUNDARY TRIPLE JUNCTIONS MODELING 


Grain boundary triple j unctions (TJs), the intersecting points of GBs, are important 
microstructural elements, especially in nanocrystalline materials. TJs may contain 
disclination and act as stress concentrators and therefore as nuclei for both plastic 
deformation and microcrack generation or formation of local amorphous regions. The 
picture on the left illustrates the structure of TJ formed by S=9/3/3 GBs in diamond, 
which was observed experimentally. Simulations indicated that several models of the TJ 
with perfect matching of the SUs can exist [1], The reconstruction along the TJ core to 
maintain four-fold coordination of atoms in diamond required double periodicity along 
the tilt axis and resulted in extended stresses in the vicinity of the TJ (fig. on the left). 

[1] O.A.Shenderova, D.W.Brenner. ‘Atomistic Simulations of Structures and 
Mechanical Properties of <01 1> Tilt Grain Boundaries and Their Triple Junctions in 
Diamond’, Phys.Rev.B 60, 7053(1999). 
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MULTISCALE MODELING OF MORE COMPLICATED 
MICROSTRUCTURAL ELEMENTS 


At the present we try to extend the multiscale DSUM model to TJs in FCC metals. 
The general strategy is the same as that used for GBs. As a first step, a wide variety of 
structural models of TJs containing special GBs should be simulated as it had been done 
by Sutton and Vitek for the structural units model for GBs. After calculations of TJs 
energies, rules for combination of SUs within a TJ core leading to a minimum energy 
core structure can be elaborated. Representing again minority SUs as disclination 
dipoles, energies and stresses of three intersecting disclination dipole walls can be 
calculated using disclination theory. TJ core energy will be obtained from atomistic 
simulation. A snapshot at the bottom illustrate a geometrical construction of a periodic 
unit cell for an idealistic nanocrystal containing only special GBs and TJs [1], Such 
system can help to understand the influence of a grain size on stability of nanocrystals. 
Again, multiscale DSUM can be applied to dramatically decrease computational burden 
of the analysis. 

[1] A.Nazarov, private communication 
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COMPUTATIONAL GRAIN BOUNDARY ENGINEERING 


The developed multiscale approach can be useful in grain boundary engineering, the 
concept that had been already applied to optimize microstructures of real materials. It is 
well known that different types of GBs could behave differently under mechanical load 
and have different resistance to crack propagation. Grain boundary engineering is aimed 
on development of a set of thermal-mechanical treatments for a polycrystalline material 
resulting in a microstructure where ‘good’ GBs (mostly low-sigma GBs) will have certain 
arrangement and proportion in a sample. Different continuum level approaches (such as 
finite elements approach) or Monte Carlo technique are used to simulate the sintering of 
polycrystalline materials and can be developed to simulate microstructure evolution 
under thermal-mechanical treatments. To make these simulations more realistic, DSUM 
could be incorporated for quick estimation of the energies of the evolving GBs and TJs in 
the simulated systems. 
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MULTISCALE MODELING OF CARBON POLYMER NANOCOMPOSITES 


A radical alternative to conventional polymer composites which has emerged during 
the last decade is a class of polymer nanocomposites (PNCs)containing a uniform 
dispersion of nanoelements with characteristic sizes less then a hundred nanometers. A 
majority of polymer molecules reside near nanoelements so that the entire matrix may be 
considered to be nanoscopically confined. These restrictions on chain conformations 
alter molecular mobility and thermal transitions. Thus the crucial issues in PNC research 
is an understanding of characteristics of the interface region, its dependence on 
nanoelement geometry and surface chemistry, the relative arrangement of constituents 
and, ultimately, its relationship to macroscopic PNC properties. While first two issues 
can be addressed by large scale atomistic simulations, the issue of collective dynamics of 
nanoelements in a matrix can be studied only using mesoscopic simulation. The last issue 
can be addresses by a traditional continuum description using characteristics provided at 
previous structural level. All level of simulations would benefit from supporting related 
experiments. 
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MESOSCOPIC DYNAMIC MODEL FOR NANOCOMPOSITE 


Wlnle large-scale atomistic studies of the interface properties had been performed 
earlier and a continuum-beam model of a nanotube was development by a group of 
researches at NASA Langley Research Center (by Harik et.al.), the model of dynamic 
behavior of an ensemble of carbon nanotubes (CNTs) in an organic matrix had not been 
considered before and is the key of our multiscale approach. In the proposed model [1], 
each nanocomposite constituent (nano-fibers and matrix) is represented through coarser- 
grained models. In particular, each CNT is modeled as a “breathing flexible cylinder” 
represented by a variable number of segments (Figure below). Segment length varies 
along the nanotube depending on the local transverse curvature. The internal interactions 
within the nanotube are described through a mesoscopic many-body force field consisting 
of terms for stretching, bending, etc. An organic matrix will be represented by a coarse- 
grained “bead-and-spring” model. The functional forms of the matrix-nanotube potentials 
are chosen based on the results of atomic-level simulations. The equations of motion for 
five independent variables are integrated and classical trajectories are obtained in a 
manner similar to the traditional MD technique. The trajectories provide complete 
information on the dynamics of the nanotubes and the matrix molecules at the 
mesoscopic length-scale. 

[1] L. Zhigilei, not published 
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PERSPECTIVES 



Multiscale Model of 
Polymer Nanocomposites: Perspectives 


The resulting general formalism for the description of CNT-polymer 
nanocomposites can be extended to other polymer nanocomposite 


materials 


Once developed and tested, the model will be used to predict the behavior and 
properties of CNT-polymer nanocomposites under conditions representative of a range of 
practical applications. In particular, mechanisms of failure during dynamic and static 
loading, material response to the shock loading, as well as acoustic and thermal 
properties of the nanocomposites will be investigated. A range of processes involved in 
fabrication and processing of nanocomposites will also be studied with the goal of 
providing recommendations for the development and optimization of fabrication and 
processing technology. 

It should be also emphasized that the general formalism for the description of 
nanocomposites implemented in the proposed model can be, in principal, extended to 
other polymer nanocomposite materials. 
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OUTLINE OF PRESENTATION 


The outline of the presentation is presented here. This includes the titles of the 
different sections of the presentation. 



Outline of Presentation 
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Introduction and Background 
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Micro-Tensile Experiments 
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Micro-Bending Experiments 


• 
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Framework for multi-scale modeling 
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Conclusions and Future Plans 
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MEMS - A MOTIVATION FOR A MULTI-SCALE FRAMEWORK 


This viewgraph identifies the emerging applications of MEMS structures. These 
include micro-motors, micro-satellites, micro-switches, accelerometers, actuators and 
sensors. The existing $ 2 billion market for MEMS structures is expected to grow to ~$ 5 
billion within a year. Photographs of emerging MEMS products are presented. These 
include an accelerometers and a turbine rotor controlled by a fiber optic cable. 


MEMS - A Motivation for a Multi-Scale 
Framework 

• Micro-electro-mechanical systems (MEMS) have emerged in 
recent years 

• These are machines on the micron-scale (finer than human hair 
of ~ 1 00 pm in diameter) 

• Emerging systems include micro-satellites, micro-motors, micro- 
switches, accelerometers, actuators and sensors 


Accelerometers Turbine Rotor 
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MATERIALS THAT ARE USED IN MEMS STRUCTURES 


This viewgraph summarizes the materials that are used currently in MEMS structures. 
Most structures are fabricated from Si (single crystal or polycrystalline). SiC is also used 
in harsh environments or high temperature environments. Most recently, metallic MEMS 
structures have emerged for larger and thicker devices e.g. Ni, Ni-Fe, Ni-Co and A1 
structures. 


Materials That Are Used in MEMS Structures 

• Most of the existing MEMS structures are 
fabricated from single crystal or poly crystal line Si 

• SiC is also used for MEMS in high 
temperature/harsh environments 

• Metallic MEMS structures are being used 
increasingly in larger and thicker devices 
(plasticity becomes important) 

- Ni, Ni-Fe, Ni-Co 

- A1 
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INTRODUCTION TO LIGA PROCESS 


This viewgraph describes the German acronym for the LIGA process that is used for 
the fabrication of the nickel MEMS structures that were employed in this study. The 
method was first developed in Germany about 20 years ago, and is now widely used for 
the fabrication of metallic MEMS structures. 


Introduction to LIGA Process 

• LIGA, an acronym of the German words for Lithogrphic, 
Galvanofonnung, Abformun (lithography, electroplating, 
and molding), is a technique used to produce micro- 
elecron-mechanical systems (MEMS) made from metals, 
ceramics, polymers and plastics 

• First developed at the Forschungszentrum Karlsruhe 
(FZK), Germany in the mid-1980s 

• Typically has 1-10 pm minimum lateral feature size, and 
are from hundreds of microns to a few millimeters thick. 
Especially advantageous for micromechanical systems 
such as micro-motors and micro-pumps because the high 
aspect ratios allow high torque to be produced 


Figure 5 
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INTRODUCTION TO LIGA PROCESS (CONT.) 


Further details on the LIGA processing of metallic MEMS structures are provided. 
The steps include the deep etching of x-ray resist materials to produce micro-moulds with 
high aspect ratios. Metals can then be electrodeposited into the micro-moulds to produce 
movable micro-devices with typical dimensions between ~ 20 - 200 pm. 


Introduction to LIGA Process (Cont.) 

• An additive microfabrication process in which structural 
material is deposited into a polymethylmethacrylate 
(PMMA) molds realized by deep x-ray photolithograpy 
(DXRL) 

Schematic of LIGA Process 
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PLASTICITY AT THE MICRON AND SUB-MICRON SCALES 


The micro-scale nature of metals MEMS structures requires a different framework for 
the modeling of plasticity due to contacts, fatigue and a wide range of loading conditions. 
Conventional J2 plasticity theory breaks down due to the size effects that become 
important at the micron scale. Extended J2 plasticity is needed along with discrete 
dislocations and atomistic approaches. 


Plasticity at The Micron and 
Sub-Micron Scales 

• There are several scenarios that require the 
constitutive modeling of plasticity at the micron- 
and sub-micron scales 

• Conventional J2 deformation theory breaks down 
in this domain since it does not contain a size- 
scale 


• Discrete dislocation and atomistic approaches 
may be needed at the smallest size scales 



Figure 7 
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PLASTICITY AT MICRON SCALE 


Metals and their alloys exhibit significant size effects at the micron scales. This is 
particularly important in the case of LIGA NI MEMS structures which have sizes 
between ~ 20 and 200 pm. Prior evidence of size effects are also presented. These 
include: evidence of micro-indentation size effects on hardness data obtained for W 
crystals, and strong evidence of size effects obtained from torsional test on copper wire. 


Plasticity at Micron Scale 

• Metals and their alloys display a strong size-dependence to 
tens of microns 

• The effect is particularly important in the characterization 
of plasticity in LIGA Ni MEMS Structures 

• Plasticity cannot be characterized with conventional 
plasticity theories that include no material length scale 



Figure 8 
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EXPERIMENTS OF PLASTICITY AT MICRON SCALE 


In general, the microhardness should be independent of the size-scale of the 
indentation. However, in the range between -1-10 pm, the hardness may become 
strongly dependent of the indentation size. Evidence of such size dependence has been 
obtained by a number of researchers. This is attributed largely to strain gradient plasticity 
phenomena that give rise to size dependent phenomena. 


Experiments of Plasticity at Micron Scale 

• Microhardness should be size-independent in the absence of a 
plasticity length scale 

• The strong size dependence of hardness at the micron-scale is 
one of the most compelling pieces of evidence for the 
existence of a plasticity length-scale 

• Strong evidence of hardness dependence observed in other 
metals 

-Atkinson (1995) 

-Ma and Clark (1995) 

-De Guzman et al. (1993) 

-Poole et al. (1996) 

-McElhaney et al. (1998) 
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MODELING OF PLASTICITY AT MICRON SCALE 


A number of researchers have made significant efforts to develop strain gradient 
plasticity theories. Two groups have emerged. The first uses phenomenological theories 
to model the effects of plastic strain gradients, while the second group uses so-called 
mechanism-based strain gradient theories based on the strain contributions from 
geometrically necessary and statistically stored dislocations. 


Modeling of Plasticity at Micron Scale 

• Significant effort to develop new plasticity theory containing 
constitutive length scale that can characterize size-dependent 
plastic deformations 

• Phenomenological plasticity theories proposed for the 
modeling of plastic strain gradients 

-Fleck and Hutchinson (1994, 1997) 

-Aifantis (1984, 1992) 

-Acharya and Bassani (1996) 

-De Borst and Muhlhaus (1992) 

• Mechanism-based strain gradient plasticity (MSG) theory 
developed based on geometrically necessary dislocation 
-Gao, Huang and Nix (1999) 

-Huang, Gao and Nix (2000) 

-Nix and Gao (1998) 
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A FRAME WORK FOR MULTI-SCALE MODELING 



Gradient Plasticity 


Continuum 


In general, a multi-scale framework is needed for the modeling of plasticity. This is 
illustrated here for a simple nano-indentation modeling framework. At the smallest 
scales, electron transitions, atomistic simulations and dislocation dynamics are needed to 
model possible pile-up and stiction phenomena. Between the micron and sub-micron 
scales, strain gradient plasticity theories are needed. Finally, continuum theories are 
useful in the regime beyond the micron scale. 


Figure 1 1 


A Framework For Multi-Scale Modeling 


At onus tics and 

Dislocation 

Dynamics 


Electron 
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MATERIALS AND PROCESSING 


The compositions of the sulfamate baths that were used in the LIGA processing of the 
electrodeposited Ni structures are presented along with the key plating conditions. Plating 
was done from a sulfamate bath to electrodeposit Ni into moulds at a current density of 
50 mA/cm 2 


Materials and Processing 

• Material used in present study was obtained from Sandia 
National Lab 

• Plating from a sulfamate bath to electro-deposit Ni in molds at 
current density of 50 mA/cm 2 


Composition and Operating Conditions 
of Nickel Sulfamate Plating Baths 


Ni((NH 2 S0 3 ) 2 «4H 2 0 

440.1 g/1 

Boric Acid 

48 g/1 

Wetting Agent 

0.2 % /vol 

Temperature 

50°C 

PH 

3. 8-4.0 


i i 


University j 


Figure 12 


410 





MICROSTRUCTURES OF UIGA NI THIN FIUM 


Cross-sectional micrographs of the complex material microstructure are presented. 
These show a refined top layer on a columnar structure with micro-scale structure. Within 
the columnar structure, the cross-sectional view from the focused ion beam image reveals 
a refined nano-scale/sub-micron structure. 


Microstructure of LIGA Ni Thin Film 


Cross-sectional View Too View (SEM Image) 



Cross-sectional View (Ion Beam Image) 



Figure 13 
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MICRO-TEXTURE OF LIGA NI THIN FILM 


The micro-texture of the as-plated LIGA nickel material is presented. The orientation 
imaging microscopy (OIM) image shows strong [001] texture. However, some “black” 
areas are also present, for which the micro-texture has not been determined. 


Micro-texture of LIGA Ni Thin Film 

Cross-sectional View 
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SCHEMATICS OF LIGA NI SPECIMENS 


Schematics of the specimens that were used in the experiments are presented. These 
include dog-bone micro-tensile specimens and rectangular specimens that were used for 
micro-bending and nano-indentation testing. 


Schematics of LIGA Ni Specimens 




Figure 15 
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MICRO-TENSILE TESTING SYSTEM 


A schematic of the micro-tensile testing system is presented. This shows the loading 
configuration, laser interferometry system for gauge strain measurement, load cell and 
uni-slide drive for load application. These all rest on an anti-vibration table. 


Micro-Tensile Testing System 

• Strain measurement: 

-Laser Interferometry 

-Optical Microscopy plus Image Analysis 

• Stress measurement 
-Load Cell 



Figure 16 
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STRESS-STRAIN CURVES OF LIGA NI THIN FILM 


The stress-strain curves obtained for LIGA Ni specimens with thickness of 50, 100 
and 200 pm are presented. These show an initial elastic regime that is followed by a 
plastic regime with limited hardening. All the three materials exhibit similar stress-strain 
behavior, i.e. no significant thickness effect in tension. 



Figure 17 
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SUMMARY OF MICRO-TENSILE PROPERTIES 


The micro-tensile properties obtained from the micro-tensile tests are presented. The 
measured differences obtained from the LIGA Ni MEMS structures with different 
thicknesses are within the limits of experimental error. Hence, there does not appear to be 
a significant effect of thickness on tensile properties. 


Summary of Micro-Tensile Properties 


Thickness 
| of Structure 
(pirn ) 

0.2% Offset 
Yield Stress 
(MPa) 

Ultimate 
Tensile 
Stress (MPa) 

Plastic 
Elongati 
on to 
Failure 

50 pirn 

405 

497 

6.6% 

1 00 pm 

475 

587 

6.2% 

200 pm | 450 547 

8.5% 
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ANALYSIS OF MICRO-TENSILE TESTS 


The absence of a size effect in micro-tension is associated with a statistical 
distribution of local short range strain gradient effects surrounding individual grains. 
These do not have a significant long range effect due to the different orientation of the 
grain boundaries. In contrast, long range applied strain gradients in bending have a 
significant size effect. 


Analysis of Micro-Tensile Tests 




No apparent effect of thickness on the tensile properties for 
LIGA Ni MEMS structure has been found. This is due to 
the lack of contribution from the st ain gradient plasticity 
phenomena 


Schematic of local 
Strain gradient 



Schematic of global 
Strain gradient 
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SCHEMATIC ILLUSTRATION OF GEOMETRICALLY NECESSARY 
DISLOCATIONS CONCEPT 


Schematic illustration of specimen and dislocation configurations are presented for 
conditions before and after bending. The schematics show that additional columns of 
atoms are needed for geometric compatibility when strain gradients are applied under 
bending loads. 
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COMPARISON OF CURRENT AND PREVIOUS TENSIUE PROPERTIES 


A comparison of the current tensile properties with prior tensile results is presented. 
The results suggest that the current yield stresses are comparable to prior results. The 
Young’s moduli are also comparable to prior results. 


Comparison of Current and Previous Tensile 

Properties 


Source 

Young’s 

Modulus 

(GPa) 

0.2% Offset 
Yield Stress 
(MPa) 

Ultimate 

Tensile 

Strength 

(MPa) 

Mazza et al. 

202 

405 

782 

Sharpe et al. 

176 

323 

555 

Hemker et al. 

175 

400 

540 

Current Study 

189 

437 

544 

Bulk Ni 

207 

59 

317 
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SCHEMATIC OF MICRO-BEND EXPERIMENT 


A schematic illustration of the micro-bend test is presented along with the 
corresponding normal moments versus strain relationships. The stages of the micro-bend 
tests are illustrated with numbers 1-3, while the corresponding moment/strain 
relationships are labeled on the schematic. 
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PROCEDURE OF MICRO-BEND EXPERIMENT 


Actual photographs of the micro-bending test are presented (before and after spring 
back). The test procedure is also described in same detail. This includes descriptions of 
Ni foil thinning with TEM preparation techniques, micro-bend testing around W fibers, 
and curvature/strain measurements with micro vision system. 


Procedure of Micro-Bend Experiment 


Ni MEMS specimens ground with TEM preparation 
techniques 

Experiments performed on specimens with thicknesses of 
25, 50, 100 and 200mm 

Micro-bending experiments performed around W fibers 
with Teflon spacers 

Elastic spring-back measured by analysis of micro-vision 
images 


Initial Configuration 


After Spring-Back 
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RESULTS FROM MICRO-BEND EXPERIMENTS 


The results from the micro-bending tests are presented. These show a significant size 
effect on plots of normalized bending moment versus surface strain. The smallest 
thicknesses (25 pm) have normalized bending moments that are 3 times those of the 
thickest (200 pm) specimens. 


Results from Micro-Bend Experiments 


Plots of normalized bending moment versus surface strain 
for all film thicknesses 

Strong evidence of size dependent effect is present in the 
micro-bending experiment 



Z 0 0.02 0.04 0.06 0.08 0.1 0.12 


Maximum strain, e b 
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SGP THEORY OF THE STRETCH AND ROTATION LENGTH SCALES 


An introduction to the strain gradient plasticity theory is presented. This includes a 
description of the strain tensor and second gradient of displacement, along with an 
expression for the effective strain in the extended J 2 deformation theory. 


SGP Theory of The Stretch and 
Rotation Length Scales 


The formulation is for a small strain and non-linear elastic 
solid, where both strain and strain gradients contribute to the 
strain energy density 

The strain tensor and the second gradient of the displacement 
(strain gradient) is defined as following: 

£ V = 2 ( u u + u j, i ) and n,k = u k,ij = SjKi + ~ 

Define the effective strain as a function of the deviatoric parts 
of the strain and strain gradient tensors as extension of 
traditional J 2 theory 

The deviatoric strain gradient tensor could be decomposed into 
3 unique, mutually orthogonal third order deviatoric tensor 
The effective strain of the deformation theory can be taken as: 


e= ~44 + WSkV'S + 44 2 ) 4 * 2) + 


m University 
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SGP THEORY OF THE STRETCH AND ROTATION LENGTH SCALES 

(CONT.) 


The complete version of the extended J 2 theory is presented initially. However since 
the XijXji is unessential in most problems, it is eliminated from consideration in the simple 
version that includes the stretch and rotational gradients. This is shown at the bottom the 
page. 


SGP Theory of The Stretch and 
Rotation Length Scales (Cont.) 


• The complete version of the J 2 SGP expression for the effective 
strain: 





Where r| i j k = |i k ij is the strain gradient generated and r|';j k is its 
deviator 

-%ij Xji is unessential in most problems (eliminated from 
consideration) 

-The effective strain is thus given by 


£ ~ ((%44) + ( l SG r >’,jk r l’ IJ k) +^ > l RGXij%fj) j 

Where -1 SG is the stretch length scale 
-1 RG is the rotation length scale 



Figure 26 
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PLASTICITY LENGTH SCALE ANALYSIS FROM MICRO-BEND 

EXPERIMENT 


The bilinear hardening law that was used to approximate the deformation behavior of 
the films is presented. A typical plot of true stress versus true strain is also presented for 
one of the LIGA foils. Typical values of the fitting parameters in the hardening 
expression are also presented. 


Plasticity Length Scale Analysis from Micro- 
bend Experiment 

• Plastic response of the foils was approximated to be rigid 
plastic with linear hardening and fit to the following relation: 

(7 = 2 0 +S p ,Ep 



Figure 27 
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PLASTICITY LENGTH SCALE ANALYSIS FROM 
MICRO-BEND EXPERIMENT (CONT.) 


An outline of the plasticity length scale analysis is presented. This includes 
expressions for the strain energy density in the presence of strain gradient, and an 
expression for the bending moment as a function of surface strain. 


Plasticity Length Scale Analysis from Micro- 
bend Experiment (Cont.) 

• For a rigid-plastic material with linear hardening, the 
strain energy density in the presence of strain gradients 
can be expressed as: 

W = |(£ f£ + 2S„) 

• The energy density is integrated through the thickness to 
obtain the total energy of the beam and differentiated 
with respect to the neutral axis curvature in order to 
express the bending moment, M, as a function of the 
surface strain, Sb, 


4 M 

EX 


2 

s. 


■yj2/3 2 +1+2J3 2 In 


• y / 2 / 7 2 +1 +1 

V2/F+1-1 


+ 


£ E 


9 3 


where fj= IJh 
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PLASTICITY LENGTH SCALE ANALYSIS FROM 
MICRO-BEND EXPERIMENT (CONT.) 


The extraction of the composite length scale parameter, l c , is discussed. The value of 
1 0 obtained for LIGA Ni foils in 5.6 pm. this compares very well with values of 1 0 of 5.3 ± 
0.2 pm obtained by Stolken and Evans. 


Plasticity Length Scale Analysis from Micro- 
bend Experiment ( Cont . ) 

• The micro-bend experiments provide a measure of l c 
which is given by: 



• The non-linear data fitting is performed using the built-in 
functions of the symbolic computation software 

MATHEMATICA 

• In current study, l c ~ 5.6 pm in LIGA Ni MEMS structures 

• This compares very well with the results of Stolken and 
Evans obtained from Ni foils for which lc =5.3±0.2 pm and 
Irg ~ 5 pm 

University 
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SET-UP OF NANO-INDENTATION EXPERIMENT 


The set-up for nano-indentation testing is presented. This includes a TnboScope 
nano-indentation system produced by Hysitron, Inc, Minneapolis, MN. The system is 
incorporated into a DI Dimension 3100 atomic force microscope (AFM) frame. Contact 
mode AFM is used to image the surfaces before and after indentation. 


Set-up of Nano-Indentation Experiment 


• The hardness measurements were performed using a 
TriboScope R (TriboScope is a registered trademark of Hysitron, 
Inc., Minneapolis, MN) Nanomechanical Testing System 
incorporated with DI Dimension 3100 AFM frame 

• Contact mode Atomic Force Microscopy scans can be obtained 
before and after the indentation tests 

• A load-displacement curve was recorded during each test 


The indenter head AFM Frame 
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PROCEDURE OF NANO- INDENTATION EXPERIMENTS 


The procedure for the nano-indentation experiments is described. This utilized North 
Star (three-sided cube comer) and Berkovich (three-sided pyramid type) indenter tips. 
Deformation was applied at a constant loading rate of 100 pN/sec, and a holding period 
of 15 seconds was applied. Consistent load-displacement plots were obtained. 


Procedure of Nano-indentation Experiments 

• Two types of indenter tips were used in the study 

- North Star indenter (a three-sided cube-comer type tip/sharp 
gradient) 

- Berkovich indenter (a three-sided pyramid type tip/weak gradient) 

• Peak load ranges between 200 pN to 11000 pN were explored in the 
study- these show very good repeatability 

• The loading rate was kept at 1000 pN/sec, and holding period of 1 5 


sec was applied 14000 
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NANO-INDENTATION SIZE DEPENDENCE 


The results from the nano-indentation tests are presented. The blunt Berkovich tip 
(small gradients) do not show a significant size effect. However, the sharp cube comer tip 
(large gradients) exhibits a strong size effect, which is associated with an increase in 
hardness by a factor of 3. 


Nano-indentation Size Dependence 


• Nano-indentation size dependence occurs due to plastic 
stretch gradients (For Shaip cube comer tip) 

• Typical stretch gradient length scale parameter measured to 
be ~ 0.25 - 1 pm in tungsten, copper and silver single 
crystal and copper poly-crystal (Begley & Hutchinson) 

• However no prior measurements for LIGA Ni Structures 
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MSG THEORY OF PLASTICITY LENGTH SCALE 


The mechanism-based strain gradient (MSG) theory is described qualitatively. The 
circular dislocation loops formed under the indenter are described. These are assumed to 
be circular loops of geometrically necessary dislocation. 


MSG Theory of Plasticity Length Scale 

• A modeling approach based on the geometrically necessary 
dislocations 

• Assumption has been made that the indentation made by rigid 
cone is accommodated by circular loops of geometrically 
necessary dislocations 




Figure 33 
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MSG THEORY OF PLASTICITY LENGTH SCALE (CONT.) 


The Taylor’s relation for the shear stress is written in terms of the separate 
contributions from geometrically necessary and statistically stored dislocations. Tabor’s 
relation and the Von Mises yield condition are also used to obtain an expression for the 
hardness. 


MSG Theory of Plasticity Length Scale (Cont.) 

• Taylor’s relation can be rewrite as following: 
t = apb 

where a is a constant to be taken as ~ 1/3 for FCC metal, p is the 
shear modulus and p s is the density of the statistically stored 
dislocations and p G is the density of the geometrically necessary 
dislocations 

• From Taylor’s Relation, apply Von Mises yielding condition of 
a=V 3x and take Tabor’s factor to be 3, we can have: 

H = 3a = 3y/3a/ibJaJl+-^=H 0 Jl + ^- 
V Ps V Ps 

Where H 0 = 3V3apbVp s 


IPa+Ps =<xpbjp 1 + — 

V Ps 
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MSG THEORY OF PLASTICITY LENGTH SCALE (CONT.) 


The density of geometrically necessary dislocations is estimated by considering the 
overall length of GND loops that are induced in a hemispherical volume defined by the 
contact radius. An estimate of the density of statistically stored dislocations is also 
obtained form the plateau hardness, Ho. 


MSG Theory of Plasticity Length Scale (Cont.) 


• The overall length of the geometric necessary dislocation loops 

induced by the indenter as modeled can be estimated by: 

. fa dr nhci 

A = 2kt 

1,0 s b 

where s = ba/h is the spacing between the loops of the GND 

• Assume all of the injected GND loops remain within the 
hemispherical volume V defined by the contact radius, a: 

V = -m 3 
3 

• So, the density of geometrically necessary dislocations 
becomes: 

A 3h 3 2 „ 

p a = — = 7 = tan 6 

° V 2ba 1 2 bh 

• Meanwhile, from the definition of H 0 , we can have: 

HI 

Ps ~ 27 a 2 n 2 b 2 _ , 

rriiicetoii U ni versitv a 
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MSG THEORY OF PLASTICITY LENGTH SCALE 


The ratio of the hardness is shown to be related to the ratio, h/h*, where h is the 
indentation depth and h is a characteristic length scale. The ratio of hardnesses is 
equivalent to the ratio of the yield strengths. 


MSG Theory of Plasticity Length Scale (Cond.) 

• Thus, a simple characteristic form for the depth 
dependence of the hardness can be obtained as: 



where H is the hardness for a given depth of indentation, 

t t H«. 

h; Ho is the hardness in the limit of infinite depth and h is 
a characteristic length scale which has the form as 
following: 

/ \ 2 

Ji ' = —ba 2 tan 2 # 

2 {Ho) 

• Alternatively, this form can be rewritten as: 
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PLASTICITY LENGTH SCALE ANALYZED USING 
CUBE CORNER INDENTER TIP DATA 


A normalized plot of (H/H (l ) 2 versus 1/h is presented for the cube-comer hardness 
values. This exhibits a linear relationship expected from the theory. Also, the values of h 
and Ho extracted from the plots of experimental data compare very well with theoretical 
predictions. 


Plasticity Length Scale Analyzed Using Cube 
Corner Indenter Tip Data 

• Plot of (H/H 0 ) 2 vs. 1/h in the case of Cube corner tip 

• A self-consistency check can be done by computing the expected 



l/hQun- 1 ) 
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h* 
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M- 
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b 
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a 
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2.60 
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73 

0.25 

0.33 

350.4 
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DETAIL OF THE PARAMETER USED IN MSG THEORY 
OF PLASTICITY LENGTH SCALE 


A framework for the estimation of strain gradients is presented. The stress ratio 
expression presented earlier is also reformulated in terms of strain gradient. 


Detail of the Parameter Used in MSG Theory 
of Plasticity Length Scale 

• From geometrically consideration in the modeling, defining 
the strain gradient as: 
tan# 


• Then we can rewrite the formulation in terms of strain 
gradient defined above: 



where a is the effective flow stress in the presence of strain 
gradient and g 0 is the flow stress in the absence of strain 
gradient 
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DETERMINATION OF STRETCH GRADIENT LENGTH SCALE PARAMETER 


An expression for the microstructural length scale parameter is presented. This has a 
value of = 0.89 pm for LIGA Ni films. From this, the Fleck-Hutchinson stretch gradient 
parameter, 1 S g is estimated to be 0.36 pm. 


Determination of Stretch Gradient Length 
Scale Parameter 


Can identify the microstructural length scale parameter as: 


/ =-b\ -Z- 

2 J 

For LIGA Ni MEMS films, the resulting value is l = 0.89/tot 
T he extracted microstructural length scale parameter is 
actually related to the material length scale parameter as: 


which can be approximately considered as the stretch 
gradient parameter since the stretch gradient plays a 
dominant role in the indentation experiment 
l™= 0. 36 pm 


I I o, 1 1 



437 





DISLOCATION PILE-UP AT SMALL SCALES 


Dislocation pile-up phenomena are shown to occur around the indents obtained at 
small scales. These suggest a need for discrete dislocation models. The challenge is really 
how to interface the discrete dislocation models with anisotropic length-scale plasticity 
theories. 


Dislocation Pile-Up At Small Scales 


• Evidence of dislocation pile-up observed in small-scale 
indentation 

• This suggests a need for discrete dislocations models (static or 
dynamic) 

• Challenge is how to mesh the discrete dislocations framework 
with anisotropic length-scale theory 
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DETERMINATION OF ROTATIONAL GRADIENT LENGTH SCALE 

PARAMETER 


The separate length-scale parameters are identified for LIGA Ni MEMS structures 
examined in this study. These give 1 0 = 5.6 pm, 1 S g = 0.36 pm and Irg = 5.58 pm. These 
are incorporated into an extended h expression for constitutive modeling of plasticity. 


Determination of Rotational Gradient Length 
Scale Parameter 

• The micro-bend experiments provide a measure of the 

composite strain gradient length scale parameter l c which is 
given by: 



• In current study, l c ~ 5.6 pm in LIGA Ni MEMS structures and 
Iso is estimated from nano-indentation study to be -0.36 pm 

• The rotational gradient parameter can be estimated to be - 5.58 
pm 

• These values can be incorporated into the extended L theory 
and used for the constitutive modeling of plasticity in Ni 
MEMS structures plated under the similar conditions: 
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SUMMARY AND CONCLUDING REMARKS 


A summary of the microstructure, mechanical properties and observed size effects is 
presented. No significant size effect was observed in tension. However, significant size 
effects were observed in the nano-indentation and micro-bend experiments. 


Summary and Concluding Remarks 

• Microstructure and mechanical properties investigated in 
LIGA Nickel MEMS structures 

• As-received structures have nano-scale gram sizes (d ~ 50 - 80 
nm) and primarily [001] micro-texture 

• No significant size effect observed under micro-tension 

• Significant size effects observed under micro-bendmg and 
nano-indentaion experiments (In the case of sharp cube comer 
tip) 

• Conventional contmuum plasticity theory without a length 
scale parameter can not explam the size dependence effect 
observed m the experiment 
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SUMMARY AND CONCLUDING REMARKS (CONT.) 


A summary of the results from the mechanism-based SGP and phenomenological 
SGP theories is presented along with the measured length scale parameters. The 
constitutive expression obtained for LIGA Ni is also presented. 


Summary and Concluding Remarks (Cont.) 

• Both phenomelogical and mechanism-based strain gradient 
plasticity theory were explored in the current study 

• Composite length-scale parameter, lc, of - 5.6 pm 
established for as-plated LIGA Ni from experiments 

• The stretch gradient parameter, 1 S g is estimated to be -0.36 
pm while the rotational gradient parameter, 1 RG , for the 
LIGA Ni films in current study is determined to be - 5.58 
pm 

• and constitutive expressions obtained for LIGA Ni: 

e = (^3 «) + (°' 362 x ^', k )+ (% x5 - 582 ^F 
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SUMMARY AND CONCLUDING REMARKS (CONT.) 


The current results are discussed within the context of a multi-scale modeling 
framework for the modeling of plastic deformation at the micron- and sub-micron scales. 
The need for anisotropic extensions to the isotropic length scale theories is also 
highlighted. 


Summary and Concluding Remarks (Cont.) 

• A framework is presented for the multi-scale modeling 
of plastic deformation at the micron and sub-micron 
scales 

• Plasticity length-scale parameters are proposed for the 
modeling of size effects at the micron-scale 

• Anisotropic extensions of the existing isotropic theory 
are needed 

• The plasticity length scale approach may be used to 
model size effects down into the sub-micron regime 
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